diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 64d4faeec1..580a54aecc 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -408,6 +408,74 @@ "name": "sum_residuals", "type": "F", "info": "Sum of residuals (mm)" + }, { + "name": "atof_region", + "type": "I", + "info": "a flag to know if a track reach S1, S2, or S3; is n if Sn is reached; is 0 otherwise" + }, { + "name": "atof_s1_x", + "type": "F", + "info": "x position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_y", + "type": "F", + "info": "y position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_z", + "type": "F", + "info": "z position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_path", + "type": "F", + "info": "total path from the beamline to the surface s1" + }, { + "name": "atof_s1_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s1" + }, { + "name": "atof_s2_x", + "type": "F", + "info": "x position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_y", + "type": "F", + "info": "y position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_z", + "type": "F", + "info": "z position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_path", + "type": "F", + "info": "total path from the beamline to the surface s2" + }, { + "name": "atof_s2_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s2" + }, { + "name": "atof_s3_x", + "type": "F", + "info": "x position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_y", + "type": "F", + "info": "y position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_z", + "type": "F", + "info": "z position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_path", + "type": "F", + "info": "total path from the beamline to the surface s3" + }, { + "name": "atof_s3_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s3" + }, { + "name": "atof_match", + "type": "B", + "info": "1 if there is a ATOF match, 0 if not" } ] }, { diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java index dead9dacdd..b276e4b605 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java @@ -7,6 +7,7 @@ import org.jlab.rec.ahdc.AHDCCluster.AHDCCluster; import org.jlab.rec.ahdc.DocaCluster.DocaCluster; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.alert.Track.Track; @@ -153,6 +154,39 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { bank.setFloat("chi2", row, (float) track.get_chi2()); bank.setFloat("sum_residuals", row, (float) track.get_sum_residuals()); + // track projection on ATOF surface S1 + Stepper stepper_s1 = track.get_ATOF_S1_stepper(); + if (stepper_s1 != null) { + bank.setFloat("atof_s1_x", row, (float) stepper_s1.y[0]); + bank.setFloat("atof_s1_y", row, (float) stepper_s1.y[1]); + bank.setFloat("atof_s1_z", row, (float) stepper_s1.y[2]); + bank.setFloat("atof_s1_path", row, (float) stepper_s1.sTot); + bank.setFloat("atof_s1_p", row, (float) stepper_s1.p()); + } + + // track projection on ATOF surface S2 + Stepper stepper_s2 = track.get_ATOF_S2_stepper(); + if (stepper_s2 != null) { + bank.setFloat("atof_s2_x", row, (float) stepper_s2.y[0]); + bank.setFloat("atof_s2_y", row, (float) stepper_s2.y[1]); + bank.setFloat("atof_s2_z", row, (float) stepper_s2.y[2]); + bank.setFloat("atof_s2_path", row, (float) stepper_s2.sTot); + bank.setFloat("atof_s2_p", row, (float) stepper_s2.p()); + } + + // track projection on ATOF surface S3 + Stepper stepper_s3 = track.get_ATOF_S3_stepper(); + if (stepper_s3 != null) { + bank.setFloat("atof_s3_x", row, (float) stepper_s3.y[0]); + bank.setFloat("atof_s3_y", row, (float) stepper_s3.y[1]); + bank.setFloat("atof_s3_z", row, (float) stepper_s3.y[2]); + bank.setFloat("atof_s3_path", row, (float) stepper_s3.sTot); + bank.setFloat("atof_s3_p", row, (float) stepper_s3.p()); + } + + bank.setInt("atof_region", row, track.get_ATOF_region()); + bank.setByte("atof_match", row, (byte) ((track.getATOFHits().size() > 0) ? 1 : 0)); + row++; } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java index 008e892104..3d9482b5a8 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java @@ -2,6 +2,7 @@ import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.detector.calib.utils.DatabaseConstantProvider; @@ -9,8 +10,9 @@ import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; +import org.jlab.rec.ahdc.KalmanFilter.KFHit; -public class Hit implements Comparable { +public class Hit implements Comparable, KFHit { private final int id; private final int superLayerId; @@ -33,6 +35,16 @@ public class Hit implements Comparable { private int trackId; //updated constructor with ADC + /** + * + * @param _Id is AHDC::adc row id + 1 + * @param _Super_layer super layer id + * @param _Layer layer id + * @param _Wire wire id + * @param _Doca distance from the timing information using the time2distance, matches {@link Hit.#time} + * @param _ADC raw ADC + * @param _Time calibrated time + */ public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca, double _ADC, double _Time) { this.id = _Id; this.superLayerId = _Super_layer; @@ -110,6 +122,7 @@ public Line3D getLine() { return wireLine; } + @Override public double getRadius() { return radius; } @@ -136,6 +149,7 @@ public double getY() { public double getPhi() {return phi;} + /** Get calibrated ADC */ public double getADC() {return adc;} public double getResidual() { @@ -146,18 +160,22 @@ public void setResidual(double resid) { this.residual = resid; } + /** Set calibrated ToT */ public void setToT(double _tot) { this.tot = _tot; } + /** Get calibrated ToT */ public double getToT() { return tot; } + /** Set calibrated ADC */ public void setADC(double _adc) { this.adc = _adc; } + /** Get raw ADC */ public double getRawADC() { return raw_adc; } @@ -174,11 +192,13 @@ public void setTrackId(int _trackId) { this.trackId = _trackId; } - public RealVector get_Vector() { + @Override + public RealVector getMeasurementVector() { return new ArrayRealVector(new double[]{this.doca}); } - public RealMatrix get_MeasurementNoise() { + @Override + public RealMatrix getMeasurementNoiseMatrix() { double mean_error = 0.471; // mm (no difference between adc and time) double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm double error_on_time = (0.4423*time + 13.7215)/(0.846038*time + 31.9867); // mm @@ -188,11 +208,50 @@ public RealMatrix get_MeasurementNoise() { //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } - // a signature for KalmanFilter.Hit_beam - public RealVector get_Vector_beam() { - return null; + // Projection function + @Override + public RealVector getProjectionVector(RealVector x) { + double d = this.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + return MatrixUtils.createRealVector(new double[]{d}); } + @Override + public RealVector getInnovationVector(RealVector x) { + RealVector measured = getMeasurementVector(); + RealVector predicted = getProjectionVector(x); + return measured.subtract(predicted); + } + + // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix getProjectionMatrix(RealVector x) { + + double ddocadx = partialProjectionMatrix(x, 0); + double ddocady = partialProjectionMatrix(x, 1); + double ddocadz = partialProjectionMatrix(x, 2); + double ddocadpx = partialProjectionMatrix(x, 3); + double ddocadpy = partialProjectionMatrix(x, 4); + double ddocadpz = partialProjectionMatrix(x, 5); + + return MatrixUtils.createRealMatrix(new double[][]{ + {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); + } + + private double partialProjectionMatrix(RealVector x, int i) { + double h = 1e-8;// in mm + RealVector x_plus = x.copy(); + RealVector x_minus = x.copy(); + + x_plus.setEntry(i, x_plus.getEntry(i) + h); + x_minus.setEntry(i, x_minus.getEntry(i) - h); + + double doca_plus = this.getProjectionVector(x_plus).getEntry(0); + double doca_minus = this.getProjectionVector(x_minus).getEntry(0); + + return (doca_plus - doca_minus) / (2 * h); + } + + @Override public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java deleted file mode 100644 index fd25c2e94d..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java +++ /dev/null @@ -1,52 +0,0 @@ -package org.jlab.rec.ahdc.KalmanFilter; - -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.RealVector; -import org.jlab.geom.prim.Line3D; -import org.jlab.geom.prim.Point3D; -import org.jlab.rec.ahdc.Hit.Hit; - -/** - * A weird object that just want to be considered as a Hit - * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() - * - * @author Mathieu Ouillon - * @author Éric Fuchey - * @author Felix Touchte Codjo - */ -public class Hit_beam extends Hit { - - private double x,y,z; - private double r,phi; - Line3D beamline; - - public Hit_beam(double x, double y , double z) { - super(0,0,0,0, 0, 0, -1); - this.x = x; - this.y = y; - this.z = z; - this.r = Math.hypot(x,y); - this.phi = Math.atan2(y,x); - beamline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis - } - - @Override - public RealVector get_Vector_beam() { - return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); - } - - @Override - public Line3D getLine() { - return beamline; - } - - @Override - public double distance(Point3D point3D) { - return this.beamline.distance(point3D).length(); - } - - @Override - public double getRadius() { - return r; - } -} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java new file mode 100644 index 0000000000..7cf2237048 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -0,0 +1,29 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Point3D; +/** + * An interface to unify the hits used the Kalman Filter (e.g AHDC hits, ATOF hits, beamline) + * + * @author Felix Touchte Codjo + */ +public interface KFHit { + public double distance(Point3D point3D); + public double getRadius(); + /** Return the measurement encoded in this KFHit */ + public RealVector getMeasurementVector(); + /** Return the measurement noise matrix for this this KFHit */ + public RealMatrix getMeasurementNoiseMatrix(); + /** Compute the measurement for a given state vector */ + public RealVector getProjectionVector(RealVector x); + /** Compute the Jacobian matrix of the {@link #getProjectionVector(RealVector)} with respect of the components of the state vector */ + public RealMatrix getProjectionMatrix(RealVector x); + + /** + * Compute the innovation by subtracting {@link #getMeasurementVector()} and {@link #getProjectionVector(RealVector)} + * @param x current state estimation + * @return return the innovation vector + */ + public RealVector getInnovationVector(RealVector x); +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 7263d6638f..08b468522a 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java @@ -2,7 +2,6 @@ import java.util.HashMap; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; @@ -22,24 +21,23 @@ public class KFitter { private RealVector stateEstimation; private RealMatrix errorCovariance; - public final Stepper stepper; + private Stepper stepper; private final Propagator propagator; private final HashMap materialHashMap; public double chi2 = 0; // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; - private double[] vertex_resolutions = {0.09,1e10}; // default values // dr^2 and dz^2 in mm^2 - public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator, final HashMap materialHashMap) { - this.stateEstimation = initialStateEstimate; - this.errorCovariance = initialErrorCovariance; - this.stepper = stepper; + public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Propagator propagator, final HashMap materialHashMap) { + this.stateEstimation = initialStateEstimate.copy(); + this.errorCovariance = initialErrorCovariance.copy(); + this.stepper = new Stepper(initialStateEstimate.toArray()); this.propagator = propagator; this.materialHashMap = materialHashMap; } - public void predict(Hit hit, boolean direction) throws Exception { + public void predict(KFHit hit, boolean direction) throws Exception { // Initialization stepper.initialize(direction); Stepper stepper1 = new Stepper(stepper.y); @@ -86,39 +84,17 @@ public void predict(Hit hit, boolean direction) throws Exception { errorCovariance = (transitionMatrix.multiply(errorCovariance.multiply(transitionMatrixT))).add(processNoise); } - public void correct(Hit hit) { - RealVector z; - RealMatrix measurementNoise; - RealMatrix measurementMatrix; - RealVector h; - // check if the hit is the beamline - if (hit.getRadius() < 1) { - // the diagonal elements are the squared errors in r, phi, z - measurementNoise = - new Array2DRowRealMatrix( - new double[][]{ - {vertex_resolutions[0], 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, vertex_resolutions[1]} - });//3x3 - measurementMatrix = H_beam(stateEstimation);//6x3 - h = h_beam(stateEstimation);//3x1 - z = hit.get_Vector_beam();//0! - } - // else, it is an AHDC hits - else { - measurementNoise = hit.get_MeasurementNoise();//1x1 - measurementMatrix = H(stateEstimation, hit);//6x1 - h = h(stateEstimation, hit);//1x1 - z = hit.get_Vector();//1x1 - } + public void correct(KFHit hit) { + + RealMatrix measurementNoise = hit.getMeasurementNoiseMatrix(); + RealMatrix measurementMatrix = hit.getProjectionMatrix(stateEstimation); RealMatrix measurementMatrixT = measurementMatrix.transpose(); // S = H * P(k) * H' + R RealMatrix S = measurementMatrix.multiply(errorCovariance).multiply(measurementMatrixT).add(measurementNoise); // Inn = z(k) - h(xHat(k)-) - RealVector innovation = z.subtract(h); + RealVector innovation = hit.getInnovationVector(stateEstimation); double chi2inc = innovation.dotProduct(MatrixUtils.inverse(S).operate(innovation)); chi2 += chi2inc; @@ -138,6 +114,7 @@ public void correct(Hit hit) { stepper.y = stateEstimation.toArray(); } + // specific to AHDC hits public double residual(Hit hit) { double d = hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); return hit.getDoca()-d; @@ -147,7 +124,7 @@ public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } - private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { + private RealMatrix F(KFHit hit, Stepper stepper1) throws Exception { double[] dfdx = subfunctionF(hit, stepper1, 0); double[] dfdy = subfunctionF(hit, stepper1, 1); @@ -160,7 +137,7 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { {dfdx[0], dfdy[0], dfdz[0], dfdpx[0], dfdpy[0], dfdpz[0]}, {dfdx[1], dfdy[1], dfdz[1], dfdpx[1], dfdpy[1], dfdpz[1]}, {dfdx[2], dfdy[2], dfdz[2], dfdpx[2], dfdpy[2], dfdpz[2]}, {dfdx[3], dfdy[3], dfdz[3], dfdpx[3], dfdpy[3], dfdpz[3]}, {dfdx[4], dfdy[4], dfdz[4], dfdpx[4], dfdpy[4], dfdpz[4]}, {dfdx[5], dfdy[5], dfdz[5], dfdpx[5], dfdpy[5], dfdpz[5]}}); } - double[] subfunctionF(Hit hit, Stepper stepper1, int i) throws Exception { + double[] subfunctionF(KFHit hit, Stepper stepper1, int i) throws Exception { double h = 1e-8;// in mm Stepper stepper_plus = new Stepper(stepper1.y); Stepper stepper_minus = new Stepper(stepper1.y); @@ -184,87 +161,6 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { return new double[]{dxdi, dydi, dzdi, dpxdi, dpydi, dpzdi}; } - // Measurement matrix in 1x1 dimension: minimize distance - doca - private RealVector h(RealVector x, Hit hit) { - double d = hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - return MatrixUtils.createRealVector(new double[]{d}); - } - - // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) - private RealMatrix H(RealVector x, Hit hit) { - - double ddocadx = subfunctionH(x, hit, 0); - double ddocady = subfunctionH(x, hit, 1); - double ddocadz = subfunctionH(x, hit, 2); - double ddocadpx = subfunctionH(x, hit, 3); - double ddocadpy = subfunctionH(x, hit, 4); - double ddocadpz = subfunctionH(x, hit, 5); - - return MatrixUtils.createRealMatrix(new double[][]{ - {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); - } - - double subfunctionH(RealVector x, Hit hit, int i) { - double h = 1e-8;// in mm - RealVector x_plus = x.copy(); - RealVector x_minus = x.copy(); - - x_plus.setEntry(i, x_plus.getEntry(i) + h); - x_minus.setEntry(i, x_minus.getEntry(i) - h); - - double doca_plus = h(x_plus, hit).getEntry(0); - double doca_minus = h(x_minus, hit).getEntry(0); - - return (doca_plus - doca_minus) / (2 * h); - } - - // Measurement matrix for the beamline (Hit_beam) in dimeansion 3x1 - private RealVector h_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - double zz = x.getEntry(2); - - double r = Math.hypot(xx, yy); - double phi = Math.atan2(yy, xx); - - return MatrixUtils.createRealVector(new double[]{r, phi, zz}); - } - - // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) - private RealMatrix H_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - - double drdx = (xx) / (Math.hypot(xx, yy)); - double drdy = (yy) / (Math.hypot(xx, yy)); - double drdz = 0.0; - double drdpx = 0.0; - double drdpy = 0.0; - double drdpz = 0.0; - - double dphidx = -(yy) / (xx * xx + yy * yy); - double dphidy = (xx) / (xx * xx + yy * yy); - double dphidz = 0.0; - double dphidpx = 0.0; - double dphidpy = 0.0; - double dphidpz = 0.0; - - double dzdx = 0.0; - double dzdy = 0.0; - double dzdz = 1.0; - double dzdpx = 0.0; - double dzdpy = 0.0; - double dzdpz = 0.0; - - return MatrixUtils.createRealMatrix( - new double[][]{ - {drdx, drdy, drdz, drdpx, drdpy, drdpz}, - {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, - {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} - }); - } public RealVector getStateEstimationVector() { return stateEstimation.copy(); @@ -273,10 +169,8 @@ public RealVector getStateEstimationVector() { public RealMatrix getErrorCovarianceMatrix() { return errorCovariance.copy(); } - - public void setVertexResolution(double[] res) { - vertex_resolutions[0] = res[0]; - vertex_resolutions[1] = res[1]; - } + + /** Return a copy of the stepper. It is like a snapshot of the propagation. */ + public Stepper getStepper() {return new Stepper(stepper);} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index 26dcf2ee7f..7380ccefd6 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java @@ -1,22 +1,24 @@ package org.jlab.rec.ahdc.KalmanFilter; import java.util.ArrayList; -import java.util.Collections; import java.util.HashMap; +import java.util.Random; +import java.util.logging.Level; +import java.util.logging.Logger; import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; -import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.io.base.DataBank; -import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.alert.Track.Track; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Component; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.geom.prim.Point3D; -//import org.apache.commons.math3.linear.RealMatrixFormat; /** * This is the main routine of the Kalman Filter. The fit is done by a KFitter @@ -31,149 +33,380 @@ public class KalmanFilter { public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;} - public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} + static final Logger LOGGER = Logger.getLogger(KalmanFilter.class.getName()); + + HashMap materialHashMap = MaterialMap.generateMaterials(); private PDGParticle particle; private int Niter = 40; // number of iterations for the Kalman Filter - private boolean IsVtxDefined = false; // implemented but not used yet - private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} - // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration - private double clas_alignement = -54; + private double vz_constraint = 0; + private boolean IsVtxDefined = false; + + // mm, they are the misalignement with respect to the AHDC: the are defined in ALERTEngine + private double atof_alignement = 0; - public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { + private int counter = 0; // number of utilisation of the Kalman Filter + + AlertTOFDetector ATOFdet = null; // reference to the ATOF geometry + HashMap> ATOF_hits_predicted = new HashMap<>(); // trackid vs (sector, layer, wedge) + + public void propagation(ArrayList tracks, final double magfield, boolean IsMC) { try { - double vz_constraint = 0; // to be linked to the electron vertex + counter++; // Initialization --------------------------------------------------------------------- final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; - HashMap materialHashMap = MaterialMap.generateMaterials(); - // Recover the vertex of the electron - if (event.hasBank("REC::Particle")) { - DataBank recBank = event.getBank("REC::Particle"); - int row = 0; - while ((!IsVtxDefined) && row < recBank.rows()) { - if (recBank.getInt("pid", row) == 11) { - IsVtxDefined = true; - vz_constraint = 10*recBank.getFloat("vz",row) - (IsMC ? 0 : clas_alignement); // mm - //////////////////////////////////////// - /// compute electron resolution here - /// it depends en p and theta - /// the fine tuning will be done later - /// //////////////////////////////////// - //double px = recBank.getFloat("px",row); - //double py = recBank.getFloat("py",row); - //double pz = recBank.getFloat("pz",row); - //double p = Math.sqrt(px*px+py*py+pz*pz); - //double theta = Math.acos(pz/p); - vertex_resolutions[0] = 0.09; - vertex_resolutions[1] = 64;//4 + 1e10*theta + 1e10*p; - } - row++; - } - } - + + // Loop over tracks for (Track track : tracks) { - // Initialize state vector + /// Initialize state vector double x0 = 0.0; double y0 = 0.0; - double z0 = IsVtxDefined ? vz_constraint : track.get_Z0(); + double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); double pz0 = track.get_pz(); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; - // Read list of hits + + /// Read list of hits + RadialKFHit beam_hit = track.getBeamlineHit(); ArrayList AHDC_hits = track.getHits(); - Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java + ArrayList ATOF_hits = track.getATOFHits(); - // Start propagation - Stepper stepper = new Stepper(y); + /// Initialize propagator RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B); Propagator propagator = new Propagator(RK4); - // Initialization of the Kalman Fitter - // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 - RealVector initialStateEstimate = new ArrayRealVector(stepper.y); - RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); - KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); - if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); + /// Initialization of the Kalman Fitter + RealVector initialStateEstimate = new ArrayRealVector(y); + RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); + KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, propagator, materialHashMap); // Loop over number of iterations for (int k = 0; k < Niter; k++) { - // Forward propagation + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); } - // Backward propagation (last layer to first layer) - for (int i = AHDC_hits.size() - 2; i >= 0; i--) { + + // Take into account ATOF hits + if (!AHDC_hits.isEmpty() && !ATOF_hits.isEmpty()) { + // Forward propagation in ATOF + for (KFHit hit : ATOF_hits) { + TrackFitter.predict(hit, true); + TrackFitter.correct(hit); + } + // Backward propagation in ATOF + for (int i = ATOF_hits.size()-2; i >= 0; i--){ + KFHit hit = ATOF_hits.get(i); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); + } + // Backward propagation to the last AHDC layer + { + Hit hit = AHDC_hits.getLast(); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); + } + } + + // Backward propagation (from AHDC last layer to first AHDC layer) + for (int i = AHDC_hits.size()-2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } - // Backward propagation (first layer to beamline) + + // Backward propagation (from first AHDC layer to beamline) { - Hit hit = new Hit_beam(0, 0, vz_constraint); - TrackFitter.predict(hit, false); - TrackFitter.correct(hit); + TrackFitter.predict(beam_hit, false); + TrackFitter.correct(beam_hit); } - } - - /*RealMatrixFormat format = - new RealMatrixFormat( - "[\n", "\n]", // matrix start/end - "[", "]", // row start/end - ",\n", // column separator - " ; ", // row separator - new java.text.DecimalFormat("0.0000") - ); - System.out.println("=====> Print error matrix"); - System.out.println(format.format(TrackFitter.getErrorCovarianceMatrix()));*/ + } // end loop over nb. iteration - + /// Output: position and momentum RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); + track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); - // Post fit propagation (no correction) to set the residuals - KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); + /// Post fit propagation (no correction) + KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Propagator(RK4), materialHashMap); + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); - if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits - hit.setResidual(PostFitPropagator.residual(hit)); - } + hit.setResidual(PostFitPropagator.residual(hit)); // output: residual } - - // Fill track and hit bank + + // Fill values for AHDC::kftrack // TO DO : s and p_drift have to be checked to be sure they represent what we want - double s = PostFitPropagator.stepper.sTot; - double p_drift = PostFitPropagator.stepper.p(); - int sum_adc = 0; + Stepper current_stepper = PostFitPropagator.getStepper(); + double s = current_stepper.sTot; + double p_drift = current_stepper.p(); double sum_residuals = 0; double chi2 = 0; for (Hit hit : AHDC_hits) { - sum_adc += hit.getADC(); sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0); + chi2 += Math.pow(hit.getResidual(),2)/hit.getMeasurementNoiseMatrix().getEntry(0,0); } - track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); track.set_chi2(chi2/(AHDC_hits.size()-3)); track.set_p_drift(p_drift); - track.set_dEdx(sum_adc/s); track.set_path(s); - track.set_n_hits(AHDC_hits.size()); + + /// At the end of the PostFit propagation towards the last layer of the AHDC + /// we project the track on the ATOF surface without any AI ATOF hit indication + if (ATOFdet != null) { + // Projection towards the ATOF surfaces + // R1 : radius of the lower surface of an ATOF bar + // R2 : radius of the upper surface of an ATOF bar = lower surface of an ATOF wedge + // R3 : radius of the upper surface of an ATOF wedge + // Illustration of an ATOF component : p0 and p3 are the points on the lower surface; p1 and p2 are the points on the upper surface + // p2 ---------- p1 + // \ / + // \ / + // p3 ------ p0 + Point3D pR1 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(0); + Point3D pR2 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(2); + Point3D pR3 = ATOFdet.getSector(0).getSuperlayer(1).getLayer(0).getComponent(0).getVolumePoint(2); + double R1 = Math.hypot(pR1.x(), pR1.y()); + double R2 = Math.hypot(pR2.x(), pR2.y()); + double R3 = Math.hypot(pR3.x(), pR3.y()); + + /// From last AHDC hit to surface R1 + RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); + PostFitPropagator.predict(hitR1, true); + Stepper stepperR1 = PostFitPropagator.getStepper(); + track.set_ATOF_S1_stepper(stepperR1); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { + track.set_ATOF_region(1); // indicate if we reach this surface or not + } + + /// From surface R1 to surface R2 + RadialSurfaceKFHit hitR2 = new RadialSurfaceKFHit(R2); + PostFitPropagator.predict(hitR2, true); + Stepper stepperR2 = PostFitPropagator.getStepper(); + track.set_ATOF_S2_stepper(stepperR2); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { + track.set_ATOF_region(2); // indicate if we reach this surface or not + } + // predict the wedge + double[] vecR2 = PostFitPropagator.getStateEstimationVector().toArray(); + Point3D posR2 = new Point3D(vecR2[0], vecR2[1], vecR2[2]); + posR2.translateXYZ(0, 0, atof_alignement); + ATOF_hits_predicted.put(track.get_trackId(), predict_adjacent_wedges(ATOFdet, posR2)); + + + /// From surface R2 to surface R3 + RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); + PostFitPropagator.predict(hitR3, true); + Stepper stepperR3 = PostFitPropagator.getStepper(); + track.set_ATOF_S3_stepper(stepperR3); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { + track.set_ATOF_region(3); // indicate if we reach this surface or not + } + + } // end propagation towards ATOF surface + }//end of loop on track candidates } catch (Exception e) { - //e.printStackTrace(); - //System.out.println("======> Kalman Filter Error"); + LOGGER.log(Level.FINE, "Kalman Filter propagation failed...", e); + } + } + public void set_Niter(int Niter) {this.Niter = Niter;} + public int get_Niter() {return this.Niter;} + public void set_particle(PDGParticle particle) {this.particle = particle;} + public PDGParticle get_particle() {return this.particle;} + public HashMap> get_ATOF_hits_predicted() {return this.ATOF_hits_predicted;} + + // Test + public static void main(String[] args) { + // ATOF detector + AlertTOFDetector atof = (new AlertTOFFactory()).createDetectorCLAS(new DatabaseConstantProvider()); + + int Npts = 1000; + int counter = 0; + int err = 0; + for (int i = 0; i < Npts; i++) { + Random rand = new Random(); + int true_sector = rand.nextInt(15); + int true_layer = rand.nextInt(4); + int true_wedge = rand.nextInt(10); + + Component comp = atof.getSector(true_sector).getSuperlayer(1).getLayer(true_layer).getComponent(true_wedge); + // top face + Point3D p0 = comp.getVolumePoint(0); + Point3D p1 = comp.getVolumePoint(1); + Point3D p2 = comp.getVolumePoint(2); + Point3D p3 = comp.getVolumePoint(3); + // bottom face + Point3D p4 = comp.getVolumePoint(4); + Point3D p5 = comp.getVolumePoint(5); + Point3D p6 = comp.getVolumePoint(6); + Point3D p7 = comp.getVolumePoint(7); + + // Random point int he current wedge volume + //Point3D pt = p0.lerp(p7, Math.random()); + //Point3D pt = comp.getMidpoint(); + double t0 = rand.nextDouble(1); + double t1 = rand.nextDouble(1-t0); + double t2 = rand.nextDouble(1-t0-t1); + double t3 = 1-t0-t1-t2; + //System.out.printf("t0 + t1 + t2 + t3 = %f\n", t0+t1+t2+t3); + double t4 = rand.nextDouble(1); + double t5 = rand.nextDouble(1-t4); + double t6 = rand.nextDouble(1-t4-t5); + double t7 = 1-t4-t5-t6; + //System.out.printf("t4 + t5 + t6 + t7 = %f\n", t4+t5+t6+t7); + double x_top = t0*p0.x() + t1*p1.x() + t2*p2.x() + t3*p3.x(); + double y_top = t0*p0.y() + t1*p1.y() + t2*p2.y() + t3*p3.y(); + double x_bot = t4*p4.x() + t5*p5.x() + t6*p6.x() + t7*p7.x(); + double y_bot = t4*p4.y() + t5*p5.y() + t6*p6.y() + t7*p7.y(); + Point3D pt_top = new Point3D(x_top, y_top, p0.z()); + Point3D pt_bot = new Point3D(x_bot, y_bot, p4.z()); + Point3D pt = pt_top.lerp(pt_bot, Math.random()); + //System.out.printf("distance from midpoint : %f\n", comp.getMidpoint().distance(pt)); + + // Test the algoritm + int[] res = KalmanFilter.predict_wedge(atof, pt); + + if (res[0] == true_sector && res[1] == true_layer && res[2] == true_wedge) { + counter++; + } else { + err++; + System.out.printf("%d) Initial wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, true_sector, true_layer, true_wedge); + System.out.printf("%d) Predicted wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, res[0], res[1], res[2]); + } + } + System.out.printf("Nb of testing : %d\n", Npts); + System.out.printf("Nb of success : %d (%.2f %%)\n", counter, 100.0*counter/Npts); + } + /** + * @param pt is defined in the center of the ATOF + */ + static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { + // find the wedge + int wedge = -1; + double dz = 1e10; + for (int c = 0; c < atof.getSector(0).getSuperlayer(1).getLayer(0).getNumComponents(); c++) { + Point3D midpoint = atof.getSector(0).getSuperlayer(1).getLayer(0).getComponent(c).getMidpoint(); + if (Math.abs(midpoint.z()-pt.z()) < dz) { + dz = Math.abs(midpoint.z()-pt.z()); + wedge = c; + } + } + if (wedge == -1) return null; + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(1).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(1).getLayer(l).getComponent(wedge).getMidpoint(); + if (midpoint.distance(pt) < d) { + d = midpoint.distance(pt); + sector = s; + layer = l; + } + } + } + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, wedge}; } } - void set_Niter(int Niter) {this.Niter = Niter;} - void set_particle(PDGParticle particle) {this.particle = particle;} + + /** + * @param pt is defined in the center of the ATOF + */ + static public int[] predict_bar(AlertTOFDetector atof, Point3D pt) { + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(0).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(0).getLayer(l).getComponent(10).getMidpoint(); + double distance = midpoint.vectorTo(pt).rho(); + if (distance < d) { + d = distance; + sector = s; + layer = l; + } + } + } + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, 10}; + } + } + + public ArrayList predict_adjacent_wedges(AlertTOFDetector atof, Point3D pt) { + int[] closest_wedge_id = predict_wedge(atof, pt); + if (closest_wedge_id == null) { + return new ArrayList(); // an empty list + } + + int sector = closest_wedge_id[0]; + int layer = closest_wedge_id[1]; + int wedge = closest_wedge_id[2]; + + // find adjacent layer and sector + int sector_plus = sector; + int sector_minus = sector; + int layer_plus = layer+1; + int layer_minus = layer-1; + if (layer == 0) { + sector_plus = sector; + sector_minus = Math.floorMod(sector-1, 15); + layer_plus = layer+1; + layer_minus = 3; + } + else if (layer == 3) { + sector_plus = Math.floorMod(sector+1, 15); + sector_minus = sector; + layer_plus = 0; + layer_minus = layer-1; + } + // Here are all the adjacents wedges (maximum 9) + ArrayList listOfWedges = new ArrayList<>(); + listOfWedges.add(new int[]{sector, layer, wedge}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge}); + if (wedge-1 >= 0) { + listOfWedges.add(new int[]{sector, layer, wedge-1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge-1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge-1}); + } + if (wedge+1 <= 9) { + listOfWedges.add(new int[]{sector, layer, wedge+1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge+1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge+1}); + } + + // order the list following the distance to the initial point + listOfWedges.sort((a,b) -> { + double da = atof.getSector(a[0]).getSuperlayer(1).getLayer(a[1]).getComponent(a[2]).getMidpoint().distance(pt); + double db = atof.getSector(b[0]).getSuperlayer(1).getLayer(b[1]).getComponent(b[2]).getMidpoint().distance(pt); + return Double.compare(da, db); + }); + + return listOfWedges; + } + + public void set_ATOF_detector(AlertTOFDetector atof) { this.ATOFdet = atof;} + public void set_atof_alignement(double _shift) {this.atof_alignement = _shift;} + public void set_vz_constraint(double _vz) {this.vz_constraint = _vz;} + public void set_vertex_flag(boolean _flag) {this.IsVtxDefined = _flag;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java index e85021addb..6f4f58a835 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java @@ -1,13 +1,10 @@ package org.jlab.rec.ahdc.KalmanFilter; -//import java.util.Arrays; - import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Point3D; import org.jlab.clas.tracking.kalmanfilter.Material; import java.util.HashMap; -import org.jlab.rec.ahdc.Hit.Hit; /** * Is responsible of the propagation @@ -27,7 +24,7 @@ public Propagator(RungeKutta4 rungeKutta4) { } // Propagate the stepper toward the next hit - void propagate(Stepper stepper, Hit hit, HashMap materialHashMap) { + void propagate(Stepper stepper, KFHit hit, HashMap materialHashMap) { // Do not allow more than 10000 steps (very critical cases) final int maxNbOfStep = 10000; @@ -95,9 +92,8 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } else { // the distance between the step is not so big but the distance with respect to the hit starts to increase - // so, go back to the previous step (the best we have) and stop the propagation, set the default step size + // so, go back to the previous step (the best we have) and stop the propagation stepper.copyContent(prev_stepper); - stepper.h = 0.5; break; } } @@ -161,7 +157,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } - public RealVector f(Stepper stepper, Hit hit, HashMap materialHashMap) { + public RealVector f(Stepper stepper, KFHit hit, HashMap materialHashMap) { propagate(stepper, hit, materialHashMap); return new ArrayRealVector(stepper.y); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java new file mode 100644 index 0000000000..2f2dbebc66 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -0,0 +1,128 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; + +/** + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ +public class RadialKFHit implements KFHit { + + private double r,phi,z; + Line3D line; + RealMatrix measurementNoise = null; + + public RadialKFHit(double x, double y , double z) { + this.z = z; + this.r = Math.hypot(x,y); + this.phi = Math.atan2(y,x); + line = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis + } + + @Override + public RealVector getMeasurementVector() { + return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); + } + + @Override + public RealMatrix getMeasurementNoiseMatrix() { + return measurementNoise; + } + + public void setMeasurementNoise(RealMatrix measurementNoise) { + this.measurementNoise= measurementNoise; + } + + @Override + public double distance(Point3D point3D) { + return this.line.distance(point3D).length(); + } + + @Override + public double getRadius() { + return r; + } + + // Projection function + @Override + public RealVector getProjectionVector(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + double zz = x.getEntry(2); + + double r = Math.hypot(xx, yy); + double phi = Math.atan2(yy, xx); + + return MatrixUtils.createRealVector(new double[]{r, phi, zz}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix getProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + double dphidx = -(yy) / (xx * xx + yy * yy); + double dphidy = (xx) / (xx * xx + yy * yy); + double dphidz = 0.0; + double dphidpx = 0.0; + double dphidpy = 0.0; + double dphidpz = 0.0; + + double dzdx = 0.0; + double dzdy = 0.0; + double dzdz = 1.0; + double dzdpx = 0.0; + double dzdpy = 0.0; + double dzdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz}, + {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, + {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} + }); + } + + @Override + public RealVector getInnovationVector(RealVector x) { + RealVector measured = getMeasurementVector(); + RealVector predicted = getProjectionVector(x); + + double measuredRadius = measured.getEntry(0); + double measuredPhi = measured.getEntry(1); + double measuredZ = measured.getEntry(2); + + double predictedRadius = predicted.getEntry(0); + double predictedPhi = predicted.getEntry(1); + double predictedZ = predicted.getEntry(2); + + double dr = measuredRadius - predictedRadius; + double dz = measuredZ - predictedZ; + double dphi = Math.atan2( + Math.sin(measuredPhi - predictedPhi), + Math.cos(measuredPhi - predictedPhi) + ); // this ensures, we obtain an angle between -pi and pi + + return MatrixUtils.createRealVector(new double[]{dr, dphi, dz}); + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java new file mode 100644 index 0000000000..c108b92f09 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java @@ -0,0 +1,85 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Point3D; + +/** + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ +public class RadialSurfaceKFHit implements KFHit { + + private double r; + + public RadialSurfaceKFHit(double r) { + this.r = r; + } + + @Override + public RealVector getMeasurementVector() { + return new ArrayRealVector(new double[] {this.r}); + } + + @Override + public RealMatrix getMeasurementNoiseMatrix() { + return new Array2DRowRealMatrix(new double[][]{{1e-8}}); + } + + @Override + public double distance(Point3D point3D) { + return Math.abs(this.r - Math.hypot(point3D.x(), point3D.y())); + } + + @Override + public double getRadius() { + return r; + } + + // Projection function + @Override + public RealVector getProjectionVector(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double r = Math.hypot(xx, yy); + + return MatrixUtils.createRealVector(new double[]{r}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix getProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz} + }); + } + + @Override + public RealVector getInnovationVector(RealVector x) { + RealVector measured = getMeasurementVector(); + RealVector predicted = getProjectionVector(x); + return measured.subtract(predicted); + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java index b192ca978e..0331dc0f93 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java @@ -103,8 +103,7 @@ public void doOneStep(Stepper stepper) { } private double[] f(double[] y) { - double charge = 1.0; // the charge should by given by the particle object - // for now, we assume it is a proton + double charge = particle.charge(); double pModuleInverse = 1.0 / Math.sqrt(y[3] * y[3] + y[4] * y[4] + y[5] * y[5]); double k = charge * PhysicsConstants.speedOfLight() * 10 * pModuleInverse; @@ -133,13 +132,14 @@ private double[] f(double[] y) { private void energyLoss( double[] yIn, double h, org.jlab.clas.tracking.kalmanfilter.Material material) { double mass = particle.mass() * 1000; //particle mass defined in GeV, converted to MeV + int charge = particle.charge(); h /= 10; // h defined in mm, converted to cm double mom = Math.sqrt(yIn[3] * yIn[3] + yIn[4] * yIn[4] + yIn[5] * yIn[5]); double E = Math.sqrt(mom * mom + mass * mass); //material::getEloss(double p, double m) uses GeV and cm //see common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Material.java - double dedx = material.getEloss(mom/1000, mass/1000) * 1000;//Momentum, mass input in GeV, output in GeV/cm, converted to MeV/cm + double dedx = material.getEloss(mom/1000, mass/1000, charge) * 1000;//Momentum, mass input in GeV, output in GeV/cm, converted to MeV/cm double DeltaE = dedx * h; stepper.dEdx += DeltaE; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java index e108ee3493..43b9634bbb 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java @@ -1,9 +1,13 @@ package org.jlab.rec.alert.Track; import org.jlab.rec.alert.AI.InterCluster; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.rec.ahdc.AHDCCluster.AHDCCluster; import org.jlab.rec.ahdc.HelixFit.HelixFitObject; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import java.util.ArrayList; import java.util.List; @@ -95,7 +99,6 @@ public double get_pz() { public double get_chi2() {return chi2;} public double get_sum_residuals() {return sum_residuals;} // AHDC::track - public void set_dEdx(double _dEdx) { dEdx = _dEdx;} public void set_p_drift(double _p_drift) { p_drift = _p_drift;} public void set_path(double _path) { path = _path;} public double get_dEdx() { @@ -119,9 +122,7 @@ public double get_dEdx() { public int get_trackId() { return candidate.get_trackId(); } public void set_trackId(int id) { candidate.set_trackId(id); } public int get_n_hits() { return candidate.get_n_hits(); } - public void set_n_hits(int n) { candidate.set_n_hits(n); } public int get_sum_adc() { return candidate.get_sum_adc(); } - public void set_sum_adc(int s) { candidate.set_sum_adc(s); } public ArrayList getHits() { return candidate.getHits(); } public List get_Clusters() { return candidate.get_Clusters(); } public List getInterclusters() { return candidate.getInterclusters(); } @@ -140,4 +141,67 @@ public double get_dEdx() { public int get_predicted_ATOF_sector() { return candidate.get_predicted_ATOF_sector(); } public int get_predicted_ATOF_layer() { return candidate.get_predicted_ATOF_layer(); } public int get_predicted_ATOF_wedge() { return candidate.get_predicted_ATOF_wedge(); } + + // Recent extension of the Kalman Filter + /** + *

Error covariance matrix of the measurment around the beamline

+ * + *

This quite useful when one plan to run the Kalman filter several times. This object stores the current status of the error covariance matrix at the end of the backward propagation of an iteration.

+ * + *

First 3 lines in mm^2; last 3 lines in MeV^2 (in the beamline)

+ */ + RealMatrix errorCovarianceMatrix = MatrixUtils.createRealMatrix(new double[][]{ + {50 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 50 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 900 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 100 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 100 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 900}}); + + /** + * Get the starting error covariance matrix {@link #errorCovarianceMatrix} (around the beamline) for the current KF iteration. + */ + public RealMatrix getErrorCovarianceMatrix() {return errorCovarianceMatrix;} + + private ArrayList ATOF_hits = new ArrayList<>(); + private RadialKFHit beamline_hit = null; + + public ArrayList getATOFHits() { + return this.ATOF_hits; + } + + public RadialKFHit getBeamlineHit() { + return this.beamline_hit; + } + + public void setATOFHits(ArrayList _ATOF_hits) {this.ATOF_hits = _ATOF_hits;} + + public void setBeamlineHit(RadialKFHit _beamline_hit) {this.beamline_hit = _beamline_hit;} + + /** + * Update the error covariance matrix {@link #errorCovarianceMatrix} + */ + public void setErrorCovarianceMatrix(RealMatrix errorCovarianceMatrix) {this.errorCovarianceMatrix = errorCovarianceMatrix;} + + // Position and momentum when the track crosses the ATOF surface + // S1 : lower surface of an ATOF bar + // S2 : upper surface of an ATOF bar = lower surface of an ATOF wedge + // S3 : upper surface of an ATOF wedge + Stepper ATOF_S1_stepper; + Stepper ATOF_S2_stepper; + Stepper ATOF_S3_stepper; + double ATOF_S1_radius; + double ATOF_S2_radius; + double ATOF_S3_radius; + int ATOF_region = 0; // is n if the trach reaches Sn, 0 otherwise (i.e does not reach S1) + + // Projection of the Track on the ATOF surfaces + public void set_ATOF_S1_stepper(Stepper _stepper) {this.ATOF_S1_stepper = _stepper;} + public void set_ATOF_S2_stepper(Stepper _stepper) {this.ATOF_S2_stepper = _stepper;} + public void set_ATOF_S3_stepper(Stepper _stepper) {this.ATOF_S3_stepper = _stepper;} + public void set_ATOF_region(int _n) {this.ATOF_region = _n;} + public Stepper get_ATOF_S1_stepper() {return this.ATOF_S1_stepper;} + public Stepper get_ATOF_S2_stepper() {return this.ATOF_S2_stepper;} + public Stepper get_ATOF_S3_stepper() {return this.ATOF_S3_stepper;} + public int get_ATOF_region() {return this.ATOF_region;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/TrackCandidate.java b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/TrackCandidate.java index ac92c293ec..5688c9f443 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/TrackCandidate.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/TrackCandidate.java @@ -29,8 +29,6 @@ public class TrackCandidate { private final ArrayList hits = new ArrayList<>(); private int trackId = -1; ///< id of the track - private int n_hits = 0; ///< number of hits - private int sum_adc = 0; ///< sum of adc (adc) /** Candidate specialization — defaults to AHDC-only; finders that build * AHDC+ATOF candidates (GNN) override it via {@link #setType}. */ @@ -136,8 +134,6 @@ public void set_trackId(int _trackId) { hit.setTrackId(_trackId); } } - public void set_n_hits(int _n_hits) { n_hits = _n_hits;} - public void set_sum_adc(int _sum_adc) { sum_adc = _sum_adc;} public int get_trackId() {return trackId;} public int get_n_hits() { if (hits == null) { diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index eac98f0cec..5f891a99c4 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -7,6 +7,7 @@ import java.util.HashSet; import java.util.Set; import java.util.concurrent.atomic.AtomicInteger; +import java.util.Iterator; import java.util.HashMap; import java.util.Map; @@ -16,6 +17,7 @@ import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.utils.groups.IndexedTable; import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; @@ -32,6 +34,7 @@ import org.jlab.rec.ahdc.DocaCluster.DocaClusterRefiner; import org.jlab.rec.ahdc.HelixFit.HelixFitJava; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; import org.jlab.rec.ahdc.TrackFindingMode; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.PreCluster.PreClusterFinder; @@ -47,6 +50,8 @@ import org.jlab.rec.alert.Track.AtofHitStub; import org.jlab.rec.alert.Track.Track; import org.jlab.rec.alert.Track.TrackCandidate; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import java.util.List; @@ -98,7 +103,7 @@ public class ALERTEngine extends ReconstructionEngine { */ private RecoBankWriter rbc; static final Logger LOGGER = Logger.getLogger(ALERTEngine.class.getName()); - Detector ATOF; // ALERT ATOF detector + private AlertTOFDetector ATOF; // ALERT ATOF detector private AlertDCDetector AHDC; // ALERT AHDC detector /** @@ -515,8 +520,42 @@ public boolean processDataEventUser(DataEvent event) { /// Kalmam Filter /// /////////////////////////////////////// - // read the list of tracks/hits from the banks AHDC::track and AHDC::hits + /// Pre conditions if (!event.hasBank("AHDC::track")) {return false;} + if (!event.hasBank("AHDC::hits")) {return false;} + + /// tmp: misalignement with respect to the center of the AHDC (mm) + double clas_alignement = +75; + double atof_alignement = 0; + + /// Read the electron vertex + double vz_electron = 0; + double[] vz_error2 = {0.09, 1e10, 1e10}; // mm^2, radians^2, mm^2, error on r, phi, z + boolean IsVertexDefined = false; + if (event.hasBank("REC::Particle")) { + DataBank recBank = event.getBank("REC::Particle"); + for (int row = 0; row < recBank.rows(); row++) { + if (recBank.getInt("pid", row) == 11) { + vz_electron = 10*recBank.getFloat("vz",row); // conversion in mm + IsVertexDefined = true; + + + //double px = recBank.getFloat("px",row); + //double py = recBank.getFloat("py",row); + //double pz = recBank.getFloat("pz",row); + //double p = Math.sqrt(px*px+py*py+pz*pz); + //double theta = Math.acos(pz/p); + + // set the resolutions on r and z! to be done + vz_error2[0] = 0.09; // should depend on p and theta + vz_error2[2] = 64; // should depend on p and theta + + break; // only look at the first electron + } + } + } + + /// Read the list of tracks/hits from the banks AHDC::track and AHDC::hits DataBank trackBank = event.getBank("AHDC::track"); DataBank hitBank = event.getBank("AHDC::hits"); ArrayList AHDC_tracks = new ArrayList<>(); @@ -562,24 +601,155 @@ public boolean processDataEventUser(DataEvent event) { newTrack.set_trackId(trackid); AHDC_tracks.add(newTrack); } - // intialise the Kalman Filter + + /// Associate the electron vertex (the beamline hit) to each track + boolean IsMC = event.hasBank("MC::Particle"); + double vz_constraint = vz_electron + (IsMC ? 0 : clas_alignement); // we don't have the misalignment in simulation + for (Track track : AHDC_tracks) { + RadialKFHit hit_beam = new RadialKFHit(0, 0, vz_constraint); + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {vz_error2[0], 0.0000 , 0.0000}, + {0.0000 , vz_error2[1], 0.0000}, + {0.0000 , 0.0000 , vz_error2[2]} + });//3x3; + hit_beam.setMeasurementNoise(measurementNoise); + track.setBeamlineHit(hit_beam); + } + + /// Look for ATOF wedge hits predicted by the AI + HashMap map_ATOF_hits = new HashMap<>(); + for (Pair pair : matched_ATOF_hit_id) { + int trackid = pair.getKey(); + int atofid = pair.getValue(); + if (trackid > 0 && atofid > 0) { + // recover the wedge + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getShort("id", row) == atofid) { + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + z += atof_alignement; // there is a shift between AHDC and ATOF (still don't know why) ! + RadialKFHit hit = new RadialKFHit(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.setMeasurementNoise(measurementNoise); + map_ATOF_hits.put(trackid, hit); + } + } + } + } + + /// Associate the ATOF hits to each track + for (Track track : AHDC_tracks) { + RadialKFHit hit = map_ATOF_hits.get(track.get_trackId()); + ArrayList list = new ArrayList<>(); + if (hit != null) list.add(hit); // for now, we only consider one hit in the ATOF + track.setATOFHits(list); + } + + /// Intialise the Kalman Filter double magfieldfactor = runBank.getFloat("solenoid", 0); double magfield = 50*magfieldfactor; - boolean IsMC = event.hasBank("MC::Particle"); PDGParticle proton = PDGDatabase.getParticleById(2212); - int Niter = 40; + int Niter = 25; KalmanFilter KF = new KalmanFilter(proton, Niter); - /////////////////////////////////////////////////////// - // first propagation : each AHDC_tracks will be fitted - /////////////////////////////////////////////////////// - KF.propagation(AHDC_tracks, event, magfield, IsMC); - ///////////////////////////////////////////// - // write the AHDC::kftrack bank in the event - ///////////////////////////////////////////// + //KF.set_ATOF_detector(null); + KF.set_ATOF_detector(ATOF); // Reference the ATOF geometry in the Kalman Filter + KF.set_atof_alignement(atof_alignement); + KF.set_vz_constraint(vz_constraint); + KF.set_vertex_flag(IsVertexDefined); + + /// Do a first propagation + KF.propagation(AHDC_tracks, magfield, IsMC); + + /// Look at the new ATOF hits predicted after projection of the track on the lower surface of the ATOF wedges + HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); + for (Track track : AHDC_tracks) { + int trackid = track.get_trackId(); + ArrayList possible_wedges = ATOF_hits_predicted.get(trackid); + boolean IsHitSelected = false; + if (possible_wedges != null) { + for (int[] id : possible_wedges) { + int sector = id[0]; + int layer = id[1]; + int wedge = id[2]; + // check if this hit exist in ATOF::hits + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("layer", row) == layer && bank_ATOFHits.getInt("component", row) == wedge) { + // create a RadialKFHit + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + z += atof_alignement; // there is a shift between AHDC and ATOF (still don't know why) ! + RadialKFHit hit = new RadialKFHit(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.setMeasurementNoise(measurementNoise); + + ArrayList list = new ArrayList<>(); + list.add(hit); // for now, we only consider one hit in the ATOF + track.setATOFHits(list); // update the list of the ATOF hit (i.e override AI ATOF hit if it exist) + + IsHitSelected = true; + break; + } // end id matching + } // end loop over atof hit + if (IsHitSelected) break; + } + } + } + + /// Clean AHDC bad hits + double sigma = 0.5; // mm + for (Track track : AHDC_tracks) { + ArrayList AHDC_hits = track.getHits(); + Iterator it = AHDC_hits.iterator(); + while (it.hasNext()) { + Hit hit = it.next(); + if (Math.abs(hit.getResidual()) > 3*sigma) { + it.remove(); + } + } + } + + /// Second propagation : each AHDC_tracks will be fitted + KF.set_Niter(15); + KF.propagation(AHDC_tracks, magfield, IsMC); + + /// Write the AHDC::kftrack bank in the event org.jlab.rec.ahdc.Banks.RecoBankWriter ahdc_writer = new org.jlab.rec.ahdc.Banks.RecoBankWriter(); DataBank recoKFTracksBank = ahdc_writer.fillAHDCKFTrackBank(event, AHDC_tracks); event.appendBank(recoKFTracksBank); - // update the AHDC::hits bank : fill the residuals + /// Update the AHDC::hits bank : fill the residuals event.removeBank("AHDC::hits"); ArrayList AHDC_hits = new ArrayList<>(); for (Track track : AHDC_tracks) {