Feat: Clean bad AHDC hits and Add ATOF wedges in the Kalman Filter

etc/bankdefs/hipo4/alert.json

@@ -408,6 +408,118 @@
         "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_s1_comp",
+        "type": "I",
+        "info": "sector*10000 + layer*100 + component : atof component on s1"
+      }, {
+        "name": "atof_s2_comp",
+        "type": "I",
+        "info": "sector*10000 + layer*100 + component : atof component on s1"
+      }, {
+        "name": "atof_s3_comp",
+        "type": "I",
+        "info": "sector*10000 + layer*100 + component : atof component on s1"
+      }, {
+        "name": "atof_s1_sigma_x",
+        "type": "F",
+        "info": "error in atof_s1_x"
+      }, {
+        "name": "atof_s1_sigma_y",
+        "type": "F",
+        "info": "error in atof_s1_y"
+      }, {
+        "name": "atof_s1_sigma_z",
+        "type": "F",
+        "info": "error in atof_s1_z"
+      }, {
+        "name": "atof_s2_sigma_x",
+        "type": "F",
+        "info": "error in atof_s2_x"
+      }, {
+        "name": "atof_s2_sigma_y",
+        "type": "F",
+        "info": "error in atof_s2_y"
+      }, {
+        "name": "atof_s2_sigma_z",
+        "type": "F",
+        "info": "error in atof_s2_z"
+      }, {
+        "name": "atof_s3_sigma_x",
+        "type": "F",
+        "info": "error in atof_s3_x"
+      }, {
+        "name": "atof_s3_sigma_y",
+        "type": "F",
+        "info": "error in atof_s3_y"
+      }, {
+        "name": "atof_s3_sigma_z",
+        "type": "F",
+        "info": "error in atof_s3_z"
       }
     ]
   }, {

reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java

@@ -1,12 +1,14 @@
 package org.jlab.rec.ahdc.Banks;
 
+import org.apache.commons.math3.linear.RealMatrix;
 import org.jlab.io.base.DataBank;
 import org.jlab.io.base.DataEvent;
 import org.jlab.rec.ahdc.AI.InterCluster;
 import org.jlab.rec.ahdc.AI.TrackPrediction;
 import org.jlab.rec.ahdc.Cluster.Cluster;
 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.ahdc.Track.Track;
 
@@ -153,6 +155,63 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> 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());
+
+			// tmp
+			bank.setInt("atof_s1_comp", row, track.get_ATOF_S1_component());
+			bank.setInt("atof_s2_comp", row, track.get_ATOF_S2_component());
+			bank.setInt("atof_s3_comp", row, track.get_ATOF_S3_component());
+
+			// error on atof_sn_(x,y,z)
+			RealMatrix atof_s1_errorMatrix = track.get_ATOF_S1_errorMatrix();
+			RealMatrix atof_s2_errorMatrix = track.get_ATOF_S2_errorMatrix();
+			RealMatrix atof_s3_errorMatrix = track.get_ATOF_S3_errorMatrix();
+			if (atof_s1_errorMatrix != null) {
+				bank.setFloat("atof_s1_sigma_x", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(0, 0)));
+				bank.setFloat("atof_s1_sigma_y", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(1, 1)));
+				bank.setFloat("atof_s1_sigma_z", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(2, 2)));
+			}
+			if (atof_s2_errorMatrix != null) {
+				bank.setFloat("atof_s2_sigma_x", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(0, 0)));
+				bank.setFloat("atof_s2_sigma_y", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(1, 1)));
+				bank.setFloat("atof_s2_sigma_z", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(2, 2)));
+			}
+			if (atof_s3_errorMatrix != null) {
+				bank.setFloat("atof_s3_sigma_x", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(0, 0)));
+				bank.setFloat("atof_s3_sigma_y", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(1, 1)));
+				bank.setFloat("atof_s3_sigma_z", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(2, 2)));
+			}
+
 			row++;
 		}
 

reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java

@@ -2,15 +2,17 @@
 
 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;
 import org.jlab.geom.detector.alert.AHDC.AlertDCDetector;
 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<Hit> {
+public class Hit implements Comparable<Hit>, KFHit {
 
 	private final int    id;
 	private final int    superLayerId;
@@ -106,10 +108,12 @@ public double getDoca() {
 		return doca;
 	}
 
+	@Override
     public Line3D getLine() {
         return wireLine;
     }
 
+	@Override
     public double getRadius() {
         return radius;
     }
@@ -174,11 +178,13 @@ public void setTrackId(int _trackId) {
 		this.trackId = _trackId;
 	}
 
-	public RealVector get_Vector() {
+	@Override
+	public RealVector MeasurementVector() {
 		return new ArrayRealVector(new double[]{this.doca});
 	}
 
-    public RealMatrix get_MeasurementNoise() {
+	@Override
+    public RealMatrix MeasurementNoiseMatrix() {
 		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 +194,43 @@ 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 ProjectionFunction(RealVector x) {
+		double d = this.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)
+	@Override
+	public RealMatrix ProjectionMatrix(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.ProjectionFunction(x_plus).getEntry(0);
+		double doca_minus = this.ProjectionFunction(x_minus).getEntry(0);
+
+		return (doca_plus - doca_minus) / (2 * h);
+	}
+
+	@Override
 	public double distance(Point3D point3D) {
 		return this.wireLine.distance(point3D).length();
 	}

reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java

@@ -0,0 +1,24 @@
+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.Line3D;
+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 Line3D getLine();
+	public double distance(Point3D point3D);
+	public double getRadius();
+	/** Return the measurement encoded in this KFHit */
+	public RealVector MeasurementVector();
+	/** Return the measurement noise matrix for this this KFHit */
+	public RealMatrix MeasurementNoiseMatrix();
+	/** Compute the measure for a given state vector */
+	public RealVector ProjectionFunction(RealVector x);
+	/** Compute the Jacobian matrix of the {@link #ProjectionFunction(RealVector)} with respect of the components of the state vector */
+	public RealMatrix ProjectionMatrix(RealVector x);
+}