New features for ALERT Kalman filter+optimized parameters (#445)

* Optimization of Kalman Filter:
* adjusted number of filtering iterations from 10 to 5;
* adjusted step size dx for calculation of ddoca/dx from 10^8 to 10^5;

* Added AHDC hits residuals (post-fit and pre-fit) in the output:
* residuals in the AHDC::Hits list in the alert.json file;
* filling the hits residuals in the RecoBankWriter;
* added residual and residual_prefit in ahdc/Hit/Hit.java
* added a identification flag to match ahdc/KalmanFilter/Hit.java to ahdc/Hit/Hit.java

* Successfully affected the calculated hit residual to the correct AHDC::Hit.

* Fixed and improved the calculation of the post-fit residuals:
* affecting the track parameters to the KFTrack right after the fit;
* redo a forward indicators pass without correction;

* * Cleaning the Kalman filter code:
- removed all "cylindrical coordinates" vector and measurement functions;
- renamed all preexisting vector and measurement function with their original name.
- removed many commented printouts.

* Attempt to include hit "sign" / left-right disambiguation:
 * added "virtual wires" located at the distance-of-closest-approach of the actual wire, on each side of the wire;
 * added hit sign parameter in KalmanFilter/Hit class;
 * added a new distance function to KalmanFilter/Hit class calculate the distance of a point to the correct virtual wire depending on the sign;
 * attempt to modify the "h" function to call new distance function

* Fix of a parameter modified by mistake.

* Added a second definition of BackwardIndicators in AHDC/KalmanFilter to
be able to initialize a vertex.

* Added a flag setDefinedVertex to AHDC/KalmanFilter and KFitter to define
"hit_beam" vertex.

* Reset Niter and ddoca step size parameters to 10, 1.e-8 respectively.

* Added reading of wire ADC from the AHDC HitReader, and functions to
access ADC for AHDC/Hit/Hit and AHDC/KalmanFilter/Hit.
Added filtering of two hits on same superlayer/layer based on ADC
(largest ADC is kept) and use info to determine the hit sign.

* Added an option to build the initial track with just the hits combination and preset fixed parameters in AHDCEngine.
Added a function in AHDC/KalmanFilter/Hit.java to calculate the measurement vector if we have a sign.

* Substituted call of default hit vector and hit measurement functions with hit vector and measurement functions that handle hit left/right disambiguation.

* Started to reintroduce the hit sign.

* Save state: back to status quo before revising sign.

* Added variable measurement error for hits with sign defined, with tracks on the wrong side.

* Implemented varaible measurement error for signed hits:
 * if track on right side of wire, normal error;
 * if track on wrong side of wire, inflated error;
Ensured reordering of hits by increasing phi; added exception for
"rollover" around phi = pi;

* fixed once and for all the convention for hit sign:
sign >0  if phi_expected state > phi_wire

* Tried to introduce a "pull" to the track on the correct sign of a wire by setting the measurement on the correct "virtual wire" with a larger error.

* Revert "Tried to introduce a "pull" to the track on the correct sign of a wire by setting the measurement on the correct "virtual wire" with a larger error."

This reverts commit 9bf4715.

* Fixed the convention for the "virtual wires": wire "minus" ("plus") at +deltaphi (-deltaphi)
since wire x, y position depend on -R*sin(phi), -R cos(phi)
respectively.

* Improved the functions to calculate hit vector: returns doca if sign is 0 or if sign is good.

* Added a hit distance function with goodsign as input, and H (measurement matrix) function with goodsign as an input.

* Added (commented) calls of functions with sign.

* Added a simple handle to disable reading of MC variables.

* Rerolled to fitting with no double hit.

* Harmonized simulation flag:
 - one simulation flag is declared in AHDCengine and defined as false;
 - it is now propagated into KalmanFilter.

* Added a check to read MC hits in AHDC_engine.

---------

Co-authored-by: Mathieu Ouillon <[email protected]>

Author: efuchey

etc/bankdefs/hipo4/alert.json

@@ -181,6 +181,14 @@
         "name": "Doca",
         "type": "D",
         "info": "distance od closest approch (mm)"
+      }, {
+	"name": "residual",
+	"type": "D",
+	"info": "residual (mm)"
+      }, {
+	"name": "residual_prefit",
+	"type": "D",
+	"info": "residual pre-fit (mm)"
       }
     ]
   }, {

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

@@ -24,6 +24,8 @@ public DataBank fillAHDCHitsBank(DataEvent event, ArrayList<Hit> hitList) {
 			bank.setByte("superlayer", i, (byte) hitList.get(i).getSuperLayerId());
 			bank.setInt("wire", i, hitList.get(i).getWireId());
 			bank.setDouble("Doca", i, hitList.get(i).getDoca());
+			bank.setDouble("residual", i, hitList.get(i).getResidual());
+			bank.setDouble("residual_prefit", i, hitList.get(i).getResidualPrefit());
 		}
 
 		return bank;

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

@@ -9,21 +9,28 @@ public class Hit implements Comparable<Hit> {
 	private final int    layerId;
 	private final int    wireId;
 	private final double doca;
+	private final double adc;
 
 	private double  phi;
 	private double  radius;
 	private int     nbOfWires;
 	private boolean use = false;
 	private double  x;
 	private double  y;
+	private double  residual_prefit;
+	private double  residual;
 
-	public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca) {
+        //updated constructor with ADC
+	public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca, double _ADC) {
 		this.id           = _Id;
 		this.superLayerId = _Super_layer;
 		this.layerId      = _Layer;
 		this.wireId       = _Wire;
 		this.doca         = _Doca;
+		this.adc          = _ADC;
 		wirePosition();
+		this.residual_prefit = 0.0;
+		this.residual        = 0.0;
 	}
 
 	private void wirePosition() {
@@ -130,4 +137,22 @@ public double getY() {
 	}
 
 	public double getPhi() {return phi;}
+
+	public double getADC() {return adc;}
+
+	public double getResidual() {
+		return residual;
+	}
+
+	public double getResidualPrefit() {
+		return residual_prefit;
+	}
+
+	public void setResidual(double resid) {
+		this.residual = resid;
+	}
+
+	public void setResidualPrefit(double resid) {
+		this.residual_prefit = resid;
+	}
 }

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

@@ -32,9 +32,10 @@ public void fetch_AHDCHits(DataEvent event) {
 				int    layer      = number % 10;
 				int    superlayer = (int) (number % 100) / 10;
 				int    wire       = bankDGTZ.getShort("component", i);
+				double adc        = bankDGTZ.getInt("ADC", i);
 				double doca       = bankDGTZ.getShort("ped", i) / 1000.0;
 
-				hits.add(new Hit(id, superlayer, layer, wire, doca));
+				hits.add(new Hit(id, superlayer, layer, wire, doca, adc));
 			}
 		}
 		this.set_AHDCHits(hits);
@@ -75,4 +76,4 @@ public void set_TrueAHDCHits(ArrayList<TrueHit> _TrueAHDCHits) {
 		this._TrueAHDCHits = _TrueAHDCHits;
 	}
 
-}
+}

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

@@ -19,13 +19,17 @@ public class Hit implements Comparable<Hit> {
 	private final double r;
 	private final double phi;
 	private final double doca;
-	private final double adc;
+	private double adc;
 	private final double numWires;
 	private final Line3D line3D;
-
-        // Comparison with:  common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java
-        // here, SuperLayer, Layer, Wire, start from 1
-        // in AlertDCFactory, same variables start from 1
+	private final Line3D line3D_plus;
+	private final Line3D line3D_minus;
+	private int    hitidx;
+	private int    hitsign;
+
+	// Comparison with:  common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java
+	// here, SuperLayer, Layer, Wire, start from 1
+	// in AlertDCFactory, same variables start from 1
 	public Hit(int superLayer, int layer, int wire, int numWire, double r, double doca) {
 		this.superLayer = superLayer;
 		this.layer      = layer;
@@ -34,6 +38,8 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do
 		this.doca       = doca;
 		this.numWires = numWire;
 		this.adc = 0;//placeholder
+		this.hitidx = -1;
+		this.hitsign = 0;
 
 		final double DR_layer = 4.0;//OK
 		final double round    = 360.0;//OK
@@ -104,42 +110,65 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do
 		Line3D wireLine = new Line3D(lPoint, rPoint);
 		//wireLine.show();
 		this.line3D = wireLine;
-	}
-    
-    //hit measurement vector in cylindrical coordinates: r, phi, z
-        public RealVector get_Vector() {
-	    // final double costhster = Math.cos(thster);
-	    // final double sinthster = Math.cos(thster);
-	    RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi(), 0});
-	    // Array2DRowRealMatrix stereo_rotation = new Array2DRowRealMatrix(new double[][]{{1, 0.0, 0.0}, {0, costhster, -sinthster}, {0, sinthster, costhster}});//rotation of wire: needed?
-	    return wire_meas;//.multiply(stereo_rotation);
+
+		//calculate the "virtual" left and right wires accounting for the DOCA 
+		double deltaphi = Math.asin(this.doca/R_layer);
+		double wx_plus     = -R_layer * Math.sin( alphaW_layer * (this.wire-1) - deltaphi );//OK
+		double wy_plus     = -R_layer * Math.cos( alphaW_layer * (this.wire-1) - deltaphi );//OK
+
+		double wx_plus_end = -R_layer * Math.sin( alphaW_layer * (this.wire-1) + thster * (Math.pow(-1, this.superLayer-1)) - deltaphi );//OK
+		double wy_plus_end = -R_layer * Math.cos( alphaW_layer * (this.wire-1) + thster * (Math.pow(-1, this.superLayer-1)) - deltaphi );//OK
+
+		line = new Line3D(wx_plus, wy_plus, -zl/2, wx_plus_end, wy_plus_end, zl/2);
+		lPoint = new Point3D();
+		rPoint = new Point3D();
+		lPlane.intersection(line, lPoint);
+		rPlane.intersection(line, rPoint);
+
+		wireLine = new Line3D(lPoint, rPoint);
+		this.line3D_plus = wireLine;
+
+		double wx_minus     = -R_layer * Math.sin( alphaW_layer * (this.wire-1) + deltaphi );//OK
+		double wy_minus     = -R_layer * Math.cos( alphaW_layer * (this.wire-1) + deltaphi );//OK
+
+		double wx_minus_end = -R_layer * Math.sin( alphaW_layer * (this.wire-1) + thster * (Math.pow(-1, this.superLayer-1)) + deltaphi );//OK
+		double wy_minus_end = -R_layer * Math.cos( alphaW_layer * (this.wire-1) + thster * (Math.pow(-1, this.superLayer-1)) + deltaphi );//OK
+
+		line = new Line3D(wx_minus, wy_minus, -zl/2, wx_minus_end, wy_minus_end, zl/2);
+		lPoint = new Point3D();
+		rPoint = new Point3D();
+		lPlane.intersection(line, lPoint);
+		rPlane.intersection(line, rPoint);
+		
+		wireLine = new Line3D(lPoint, rPoint);
+		this.line3D_minus = wireLine;
+		
 	}
 
-    //hit measurement vector in 1 dimension: minimize distance - doca
-        public RealVector get_Vector_simple() {
+        //hit measurement vector in 1 dimension: minimize distance - doca
+        public RealVector get_Vector() {
 		return new ArrayRealVector(new double[]{this.doca});
 	}
 
-    //hit measurement vector in 1 dimension: minimize distance - doca - adds hit "sign"
-        public RealVector get_Vector_sign(int sign) {
- 	    // Attempt: multiply doca by sign
-		return new ArrayRealVector(new double[]{sign*this.doca});
+        //hit measurement vector in 1 dimension with sign: if sign = 0, return doca, otherwise return 0
+        public RealVector get_Vector(int sign, boolean goodsign) {
+		if(sign == 0 || goodsign){
+			return new ArrayRealVector(new double[]{this.doca});
+		}else{
+			return new ArrayRealVector(new double[]{0.0});
+		}
 	}
 
-        public RealMatrix get_MeasurementNoise() {
-	    final double costhster = Math.cos(thster);
-	    final double sinthster = Math.cos(thster);
-	    //dR = 0.1m dphi = pi dz = L/2 
-	    Array2DRowRealMatrix wire_noise = new Array2DRowRealMatrix(new double[][]{{0.1, 0.0, 0.0}, {0.0, Math.atan(0.1/this.r), 0.0}, {0.0, 0.0, 150.0/costhster}});//uncertainty matrix in wire coordinates
-	    Array2DRowRealMatrix stereo_rotation = new Array2DRowRealMatrix(new double[][]{{1, 0.0, 0.0}, {0, costhster, -sinthster}, {0, sinthster, costhster}});//rotation of wire
-	    wire_noise.multiply(stereo_rotation);
-	    
-	    return wire_noise.multiply(wire_noise);
-	    //
+    	public RealMatrix get_MeasurementNoise() {
+		return new Array2DRowRealMatrix(new double[][]{{0.0225}});
 	}
-
-    	public RealMatrix get_MeasurementNoise_simple() {
-	    return new Array2DRowRealMatrix(new double[][]{{0.01}});
+    
+    	public RealMatrix get_MeasurementNoise(boolean goodsign) {
+	    if(goodsign){
+		return new Array2DRowRealMatrix(new double[][]{{0.0225}});
+	    }else{
+		return new Array2DRowRealMatrix(new double[][]{{2*this.doca*this.doca}});
+	    }
 	}
 
 	public double doca() {
@@ -151,8 +180,6 @@ public double doca() {
         public double phi()    {return phi;}//at z = 0;
 
         public double phi(double z)    {
-	    // double x_0 = r*sin(phi);
-	    // double y_0 = r*cos(phi);
 	    double x_z = r*Math.sin( phi + thster * z/(zl*0.5) * (Math.pow(-1, this.superLayer-1)) );
 	    double y_z = r*Math.cos( phi + thster * z/(zl*0.5) * (Math.pow(-1, this.superLayer-1)) );
 	    return Math.atan2(x_z, y_z);
@@ -161,10 +188,16 @@ public double phi(double z)    {
 	public Line3D line() {return line3D;}
 
 	public double distance(Point3D point3D) {
-	    //System.out.println("Calculating distance: ");
-	    //this.line3D.show();
-	    //point3D.show();
-	    //System.out.println(" d = " + this.line3D.distance(point3D).length());
+		return this.line3D.distance(point3D).length();
+	}
+
+	public double distance(Point3D point3D, int sign, boolean goodsign) {
+		//if(sign!=0)
+	    	//System.out.println(" r " + this.r +  " phi " + this.phi + " doca " + this.doca + " sign " + sign + " distance " + this.line3D.distance(point3D).length() + " (sign 0) " + this.line3D_plus.distance(point3D).length() + " (sign+) " + this.line3D_minus.distance(point3D).length() + " (sign-) ");
+		if(!goodsign){	    
+			if(sign>0)return this.line3D_plus.distance(point3D).length();
+			if(sign<0)return this.line3D_minus.distance(point3D).length();
+		}
 		return this.line3D.distance(point3D).length();
 	}
 
@@ -211,12 +244,33 @@ public double getADC() {
 		return adc;
 	}
 
+	public void setADC(double _adc) {
+		this.adc = _adc;
+	}
+
 	public Line3D getLine3D() {
 		return line3D;
 	}
 
 	public double getNumWires() {
 		return numWires;
 	}
+
+	public int getHitIdx() {
+		return hitidx;
+	}
+
+	public void setHitIdx(int idx) {
+		this.hitidx = idx;
+	}
+
+	public int getSign() {
+		return hitsign;
+	}
+
+	public void setSign(int sign) {
+		this.hitsign = sign;
+	}
+    
 }
 

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

@@ -9,7 +9,7 @@ public class Hit_beam extends Hit {
 	double r,phi;
 
 	public Hit_beam(int superLayer, int layer, int wire, int numWire, double doca, double x, double y , double z) {
-		super(0, 0, 0, 0, Math.hypot(x,y), 0);
+	    super(0, 0, 0, 0, Math.hypot(x,y), 0);
 		this.x = x;
 		this.y = y;
 		this.z = z;