From f35dcf1246522e20abddd973adad96cfe1aa661c Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Thu, 23 Jan 2025 15:23:56 -0500 Subject: [PATCH 01/28] 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; --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 52 +++++++++---------- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 18 +++++-- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 6 +-- 3 files changed, 42 insertions(+), 34 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index fa4d1dc477..f7093f29c6 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -106,37 +106,35 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do 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); - } - - //hit measurement vector in 1 dimension: minimize distance - doca + // //hit measurement vector in cylindrical coordinates: r, phi, z + // public RealVector get_Vector() { + // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi(), 0.0}); + // // 2-dim matrices + // // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi()}); + // return wire_meas; + // } + + // public RealVector get_Vector_z(double z) { + // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi(z), z}); + // return wire_meas; + // } + + //hit measurement vector in 1 dimension: minimize distance - doca public RealVector get_Vector_simple() { 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}); - } - - 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() { + // final double costhster = Math.cos(thster); + // final double sinthster = Math.cos(thster); + // //dR = 0.1m dphi = pi dz = L/2 + // // 3-dim matrices + // Array2DRowRealMatrix wire_noise = new Array2DRowRealMatrix(new double[][]{{0.01, 0.0, 0.0}, {0.0, Math.pow( (phi(-zl*0.5)+phi(zl*0.5))*0.5 , 2), 0.0}, {0.0, 0.0, zl*zl*0.25}} );//uncertainty matrix in wire coordinates + // // trying 2-dim matrices + // //Array2DRowRealMatrix wire_noise = new Array2DRowRealMatrix(new double[][]{{0.01, 0.0}, {0.0, Math.pow( (phi(-zl*0.5)+phi(zl*0.5))*0.5 , 2)} );//uncertainty matrix in wire coordinates + // return wire_noise; + // // + // } public RealMatrix get_MeasurementNoise_simple() { return new Array2DRowRealMatrix(new double[][]{{0.01}}); 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 3feddf5b27..1ed27be14c 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 @@ -94,9 +94,10 @@ public void correct(Indicator indicator) { h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { + // double wire_sg = wire_sign(indicator); measurementNoise = indicator.hit.get_MeasurementNoise_simple();//1x1 measurementMatrix = H_simple(stateEstimation, indicator);//6x1 - h = h_simple(stateEstimation, indicator);//.multiply(wire_sign_mat(indicator));//1x1 + h = h_simple(stateEstimation, indicator);//1x1 z = indicator.hit.get_Vector_simple();//1x1 // measurementNoise = indicator.hit.get_MeasurementNoise();//3x3 @@ -151,7 +152,10 @@ public double wire_sign(Indicator indicator) {//let's decide: positive when (ph // System.out.println(" phi state " + phi_state + " phi wire " + phi_wire);// + " phi state alt? " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))); // return MatrixUtils.createRealMatrix(new double[][]{{(phi_state-phi_wire)/Math.abs(phi_state-phi_wire)}}); // } - + public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ + this.errorCovariance = initialErrorCovariance; + } + private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { double[] dfdx = subfunctionF(indicator, stepper1, 0); @@ -166,7 +170,7 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { } double[] subfunctionF(Indicator indicator, Stepper stepper1, int i) throws Exception { - double h = 1e-8; + double h = 1e-5;//1e-8;// in mm Stepper stepper_plus = new Stepper(stepper1.y); Stepper stepper_minus = new Stepper(stepper1.y); @@ -195,6 +199,7 @@ private RealVector h(RealVector x, Indicator indicator) { double xx = x.getEntry(0); double yy = x.getEntry(1); return MatrixUtils.createRealVector(new double[]{Math.hypot(xx, yy), Math.atan2(yy, xx), x.getEntry(2)}); + //return MatrixUtils.createRealVector(new double[]{Math.hypot(xx, yy), Math.atan2(yy, xx)}); } //measurement matrix in 1 dimension: minimize distance - doca @@ -203,6 +208,11 @@ private RealVector h_simple(RealVector x, Indicator indicator) { return MatrixUtils.createRealVector(new double[]{d});//would need to have this 3x3 } + // private RealVector h_simple_sign(RealVector x, Indicator indicator, double sign) { + // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + // return MatrixUtils.createRealVector(new double[]{d*sign});//would need to have this 3x3 + // } + //measurement noise matrix in cylindrical coordinates: r, phi, z private RealMatrix H(RealVector x, Indicator indicator) { // dphi/dx @@ -254,7 +264,7 @@ private RealMatrix H_simple(RealVector x, Indicator indicator) { } double subfunctionH(RealVector x, Indicator indicator, int i) { - double h = 1e-8; + double h = 1e-5;//1e-8;// in mm RealVector x_plus = x.copy(); RealVector x_minus = x.copy(); 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 2cd0bc5593..3f4e1f5238 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 @@ -141,7 +141,6 @@ private void propagation(ArrayList tracks, DataEvent event) { RealVector initialStateEstimate = new ArrayRealVector(stepper.y); //first 3 lines in cm^2; last 3 lines in MeV^2 RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); - KFitter kFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); /* @@ -173,10 +172,11 @@ private void propagation(ArrayList tracks, DataEvent event) { // } // } - for (int k = 0; k < 10; k++) { + for (int k = 0; k < 5; k++) { //System.out.println("--------- ForWard propagation !! ---------"); - + //Reset error covariance: + //kFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { kFitter.predict(indicator); //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); From e71cf130d891b5c64bcd4d6dcaba94ea0bc7ec00 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 10 Feb 2025 21:28:45 -0500 Subject: [PATCH 02/28] 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 --- etc/bankdefs/hipo4/alert.json | 8 +++++ parent/pom.xml | 4 +-- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 2 ++ .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 20 +++++++++++++ .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 19 +++++++++--- .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 2 +- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 29 ++++++++++++------- .../java/org/jlab/rec/ahdc/Track/Track.java | 1 - 8 files changed, 67 insertions(+), 18 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 7c125eb9cc..1c524d73e1 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -25,6 +25,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)" } ] }, { diff --git a/parent/pom.xml b/parent/pom.xml index a02664fa3f..bf7929134c 100644 --- a/parent/pom.xml +++ b/parent/pom.xml @@ -62,7 +62,7 @@ ai.djl model-zoo - 0.30.0 + 0.31.1 compile @@ -90,7 +90,7 @@ ai.djl.pytorch pytorch-jni - 2.4.0-0.30.0 + 2.5.1-0.31.1 runtime 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 72b9b5f2fe..16e880dcb2 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 @@ -24,6 +24,8 @@ public DataBank fillAHDCHitsBank(DataEvent event, ArrayList 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; 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 6ee4bab6a0..95a202b5e6 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 @@ -16,6 +16,8 @@ public class Hit implements Comparable { 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) { this.id = _Id; @@ -24,6 +26,8 @@ public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca) { this.wireId = _Wire; this.doca = _Doca; wirePosition(); + this.residual_prefit = 0.0; + this.residual = 0.0; } private void wirePosition() { @@ -130,4 +134,20 @@ public double getY() { } public double getPhi() {return phi;} + + 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; + } } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index f7093f29c6..79e764d99d 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -22,10 +22,11 @@ public class Hit implements Comparable { private final double adc; private final double numWires; private final Line3D line3D; + private int hitidx; - // 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 + // 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 +35,7 @@ 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; final double DR_layer = 4.0;//OK final double round = 360.0;//OK @@ -54,7 +56,7 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do Vector3D n2 = new Vector3D(0, 0, 1); //n2.rotateY(thopen); //n2.rotateZ(thtilt); - Plane3D rPlane = new Plane3D(p2, n2);//OK + Plane3D rPlane = new Plane3D(p2, n2);//OK switch (this.superLayer) {//OK case 1: @@ -216,5 +218,14 @@ public Line3D getLine3D() { public double getNumWires() { return numWires; } + + public int getHitIdx() { + return hitidx; + } + + public void setHitIdx(int idx) { + this.hitidx = idx; + } + } 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 index e1160642e4..d3221b2939 100644 --- 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 @@ -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; 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 3f4e1f5238..324b7d7477 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 @@ -33,6 +33,8 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event) {propagation(tracks, event);} + private final int Niter = 5; + private void propagation(ArrayList tracks, DataEvent event) { try { @@ -96,15 +98,20 @@ private void propagation(ArrayList tracks, DataEvent event) { // EPAF: *the line below is for TEST ONLY!!!* //double[] y = new double[]{vxmc, vymc, vzmc, pxmc, pymc, pzmc}; //System.out.println("y = " + vxmc + ", " + vymc + ", " + vzmc + ", " + pxmc + ", " + pymc + ", " + pzmc + "; p = " + java.lang.Math.sqrt(pxmc*pxmc+pymc*pymc+pzmc*pzmc)); - // Initialization hit //System.out.println("tracks = " + tracks); ArrayList AHDC_hits = tracks.get(0).getHits(); ArrayList KF_hits = new ArrayList<>(); + + int hit_global_idx[] = {0, 0, 0, 0, 0, 0, 0, 0}; + double[] residuals = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + double[] residuals_prefit = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + nhits = 0; for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { //System.out.println("Superlayer = " + AHDC_hit.getSuperLayerId() + ", Layer " + AHDC_hit.getLayerId() + ", Wire " + AHDC_hit.getWireId() + ", Nwires " + AHDC_hit.getNbOfWires() + ", Radius " + AHDC_hit.getRadius() + ", DOCA " + AHDC_hit.getDoca()); - Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); - + Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); + hit.setHitIdx(AHDC_hit.getId()); + hit_global_idx[nhits++] = AHDC_hit.getId(); // Do delete hit with same radius // boolean aleardyHaveR = false; // for (Hit o: KF_hits){ @@ -166,13 +173,13 @@ private void propagation(ArrayList tracks, DataEvent event) { //Print out hit residuals *before* fit: // for (Indicator indicator : forwardIndicators) { - // kFitter.predict(indicator); - // if (indicator.haveAHit()) { - // System.out.println(" Pre-fit: indicator R " + indicator.R + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum() + " residual: " + kFitter.residual(indicator) + " sign " + kFitter.wire_sign(indicator) ); - // } + // kFitter.predict(indicator); + // if (indicator.haveAHit()) { + // //AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); + // } // } - for (int k = 0; k < 5; k++) { + for (int k = 0; k < Niter; k++) { //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: @@ -181,6 +188,7 @@ private void propagation(ArrayList tracks, DataEvent event) { kFitter.predict(indicator); //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); if (indicator.haveAHit()) { + //if(k==0 && indicator.hit.getHitIdx() )residuals[] = kFitter.residual(indicator);//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); kFitter.correct(indicator); //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); @@ -196,15 +204,16 @@ private void propagation(ArrayList tracks, DataEvent event) { //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); kFitter.correct(indicator); //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); + //if(k==Niter-1 && indicator.hit.getHitIdx()>=0)AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); } } } - // //Print out residuals *after* fit: + //Print out residuals *after* fit: // for (Indicator indicator : forwardIndicators) { // kFitter.predict(indicator); // if (indicator.haveAHit()) { - // System.out.println(" Post-fit: indicator R " + indicator.R + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum() + " residual: " + kFitter.residual(indicator) + " sign " + kFitter.wire_sign(indicator) ); + // //AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); // } // } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java index 8a968ba512..fee40d2b11 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java @@ -15,7 +15,6 @@ public class Track { private List _Clusters = new ArrayList<>(); private boolean _Used = false; private final ArrayList hits = new ArrayList<>(); - ; private double x0 = 0; private double y0 = 0; From b042b022b74bb5ca840a653a0704a2a9434c6873 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 10 Feb 2025 22:37:41 -0500 Subject: [PATCH 03/28] Successfully affected the calculated hit residual to the correct AHDC::Hit. --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 48 ++++++++----------- 1 file changed, 21 insertions(+), 27 deletions(-) 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 324b7d7477..c238a13d95 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 @@ -103,15 +103,10 @@ private void propagation(ArrayList tracks, DataEvent event) { ArrayList AHDC_hits = tracks.get(0).getHits(); ArrayList KF_hits = new ArrayList<>(); - int hit_global_idx[] = {0, 0, 0, 0, 0, 0, 0, 0}; - double[] residuals = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - double[] residuals_prefit = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - nhits = 0; for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { //System.out.println("Superlayer = " + AHDC_hit.getSuperLayerId() + ", Layer " + AHDC_hit.getLayerId() + ", Wire " + AHDC_hit.getWireId() + ", Nwires " + AHDC_hit.getNbOfWires() + ", Radius " + AHDC_hit.getRadius() + ", DOCA " + AHDC_hit.getDoca()); - Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); - hit.setHitIdx(AHDC_hit.getId()); - hit_global_idx[nhits++] = AHDC_hit.getId(); + Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); + hit.setHitIdx(AHDC_hit.getId()); // Do delete hit with same radius // boolean aleardyHaveR = false; // for (Hit o: KF_hits){ @@ -122,7 +117,7 @@ private void propagation(ArrayList tracks, DataEvent event) { // if (!aleardyHaveR) KF_hits.add(hit); } - + //AHDC_hits = tracks.get(0).getHits(); /* Writer hitsWiresWriter = new FileWriter("hits_wires.dat"); @@ -170,14 +165,6 @@ private void propagation(ArrayList tracks, DataEvent event) { } writer_back.close(); */ - - //Print out hit residuals *before* fit: - // for (Indicator indicator : forwardIndicators) { - // kFitter.predict(indicator); - // if (indicator.haveAHit()) { - // //AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); - // } - // } for (int k = 0; k < Niter; k++) { @@ -188,7 +175,15 @@ private void propagation(ArrayList tracks, DataEvent event) { kFitter.predict(indicator); //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); if (indicator.haveAHit()) { - //if(k==0 && indicator.hit.getHitIdx() )residuals[] = kFitter.residual(indicator);//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); + if(k==0 && indicator.hit.getHitIdx()>0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ + //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); + AHDC_hit.setResidualPrefit(kFitter.residual(indicator)); + } + } + //System.out.println( "AHDC hit corresponding index " + indicator.hit.getHitIdx() + " pre-fit residual " + kFitter.residual(indicator) );//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); + } //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); kFitter.correct(indicator); //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); @@ -204,18 +199,18 @@ private void propagation(ArrayList tracks, DataEvent event) { //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); kFitter.correct(indicator); //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); - //if(k==Niter-1 && indicator.hit.getHitIdx()>=0)AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); + if(k==Niter-1 && indicator.hit.getHitIdx()>=0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ + //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); + AHDC_hit.setResidual(kFitter.residual(indicator)); + } + } + //System.out.println( "AHDC hit corresponding index " + indicator.hit.getHitIdx() + " post-fit residual " + kFitter.residual(indicator) );//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); + } } } } - - //Print out residuals *after* fit: - // for (Indicator indicator : forwardIndicators) { - // kFitter.predict(indicator); - // if (indicator.haveAHit()) { - // //AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); - // } - // } /* Writer writer_last = new FileWriter("track_last.dat"); @@ -226,7 +221,6 @@ private void propagation(ArrayList tracks, DataEvent event) { writer_last.close(); */ - RealVector x_out = kFitter.getStateEstimationVector(); tracks.get(0).setPositionAndMomentumForKF(x_out); From 2d7ef041e00177bb67278c4380ebb102944d1e30 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 10 Feb 2025 23:02:41 -0500 Subject: [PATCH 04/28] 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; --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 36 +++++++++++-------- 1 file changed, 21 insertions(+), 15 deletions(-) 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 c238a13d95..4663d88394 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 @@ -168,14 +168,14 @@ private void propagation(ArrayList tracks, DataEvent event) { for (int k = 0; k < Niter; k++) { - //System.out.println("--------- ForWard propagation !! ---------"); - //Reset error covariance: - //kFitter.ResetErrorCovariance(initialErrorCovariance); + //System.out.println("--------- ForWard propagation !! ---------"); + //Reset error covariance: + //kFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { kFitter.predict(indicator); //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); if (indicator.haveAHit()) { - if(k==0 && indicator.hit.getHitIdx()>0){ + if( k==0 && indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); @@ -198,20 +198,11 @@ private void propagation(ArrayList tracks, DataEvent event) { if (indicator.haveAHit()) { //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); kFitter.correct(indicator); - //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); - if(k==Niter-1 && indicator.hit.getHitIdx()>=0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ - //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); - AHDC_hit.setResidual(kFitter.residual(indicator)); - } - } - //System.out.println( "AHDC hit corresponding index " + indicator.hit.getHitIdx() + " post-fit residual " + kFitter.residual(indicator) );//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidual(kFitter.residual(indicator)); - } + //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); } } } - + /* Writer writer_last = new FileWriter("track_last.dat"); for (Indicator indicator : forwardIndicators) { @@ -224,6 +215,21 @@ private void propagation(ArrayList tracks, DataEvent event) { RealVector x_out = kFitter.getStateEstimationVector(); tracks.get(0).setPositionAndMomentumForKF(x_out); + //Residual calcuation post fit: + for (Indicator indicator : forwardIndicators) { + kFitter.predict(indicator); + if (indicator.haveAHit()) { + if( indicator.hit.getHitIdx()>0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ + //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); + AHDC_hit.setResidual(kFitter.residual(indicator)); + } + } + } + } + } + //System.out.println("y_final = " + x_out + " p_final = " + kFitter.getMomentum()); } catch (Exception e) { // e.printStackTrace(); From 217dba706dd6a4cbb25c6973c7b89cdbdaad32d4 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Tue, 11 Feb 2025 17:08:51 -0500 Subject: [PATCH 05/28] * 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. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 37 +------ .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 99 +++---------------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 91 ++--------------- 3 files changed, 28 insertions(+), 199 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 79e764d99d..f9460b7080 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -56,7 +56,7 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do Vector3D n2 = new Vector3D(0, 0, 1); //n2.rotateY(thopen); //n2.rotateZ(thtilt); - Plane3D rPlane = new Plane3D(p2, n2);//OK + Plane3D rPlane = new Plane3D(p2, n2);//OK switch (this.superLayer) {//OK case 1: @@ -107,38 +107,13 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do //wireLine.show(); this.line3D = wireLine; } - - // //hit measurement vector in cylindrical coordinates: r, phi, z - // public RealVector get_Vector() { - // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi(), 0.0}); - // // 2-dim matrices - // // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi()}); - // return wire_meas; - // } - - // public RealVector get_Vector_z(double z) { - // RealVector wire_meas = new ArrayRealVector(new double[]{this.r(), this.phi(z), z}); - // return wire_meas; - // } //hit measurement vector in 1 dimension: minimize distance - doca - public RealVector get_Vector_simple() { + public RealVector get_Vector() { return new ArrayRealVector(new double[]{this.doca}); } - // public RealMatrix get_MeasurementNoise() { - // final double costhster = Math.cos(thster); - // final double sinthster = Math.cos(thster); - // //dR = 0.1m dphi = pi dz = L/2 - // // 3-dim matrices - // Array2DRowRealMatrix wire_noise = new Array2DRowRealMatrix(new double[][]{{0.01, 0.0, 0.0}, {0.0, Math.pow( (phi(-zl*0.5)+phi(zl*0.5))*0.5 , 2), 0.0}, {0.0, 0.0, zl*zl*0.25}} );//uncertainty matrix in wire coordinates - // // trying 2-dim matrices - // //Array2DRowRealMatrix wire_noise = new Array2DRowRealMatrix(new double[][]{{0.01, 0.0}, {0.0, Math.pow( (phi(-zl*0.5)+phi(zl*0.5))*0.5 , 2)} );//uncertainty matrix in wire coordinates - // return wire_noise; - // // - // } - - public RealMatrix get_MeasurementNoise_simple() { + public RealMatrix get_MeasurementNoise() { return new Array2DRowRealMatrix(new double[][]{{0.01}}); } @@ -151,8 +126,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 +134,6 @@ 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(); } 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 1ed27be14c..9039911ece 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 @@ -73,8 +73,6 @@ public void predict(Indicator indicator) throws Exception { } public void correct(Indicator indicator) { - //System.out.println(" state before: (" + stateEstimation.getEntry(0) + ", " + stateEstimation.getEntry(1) + ", " + stateEstimation.getEntry(2) + ", " + stateEstimation.getEntry(3) + ", " + stateEstimation.getEntry(4) + ", " + stateEstimation.getEntry(5) + ");" ); - //System.out.println(" state radius before: " + Math.sqrt( Math.pow(stateEstimation.getEntry(0), 2) + Math.pow(stateEstimation.getEntry(1), 2) ) ); RealVector z, z_plus, z_minus; RealMatrix measurementNoise; RealMatrix measurementMatrix; @@ -83,30 +81,19 @@ public void correct(Indicator indicator) { measurementNoise = new Array2DRowRealMatrix( new double[][]{ - // {9.00, 0.0000, 0.0000}, - // {0.00, 1e10, 0.0000}, - // {0.00, 0.0000, 1e10} - {0.09, 0.0000, 0.0000}, - {0.00, 1.e10, 0.0000}, - {0.00, 0.0000, 1.e10} + {9.00, 0.0000, 0.0000}, + {0.00, 1e10, 0.0000}, + {0.00, 0.0000, 1e10} });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { // double wire_sg = wire_sign(indicator); - measurementNoise = indicator.hit.get_MeasurementNoise_simple();//1x1 - measurementMatrix = H_simple(stateEstimation, indicator);//6x1 - h = h_simple(stateEstimation, indicator);//1x1 - z = indicator.hit.get_Vector_simple();//1x1 - - // measurementNoise = indicator.hit.get_MeasurementNoise();//3x3 - // measurementMatrix = H(stateEstimation, indicator);//6x3 - // h = h(stateEstimation, indicator);//3x1 - // z = indicator.hit.get_Vector();//3x1 - - //System.out.println(" h: r " + h.getEntry(0) + " phi " + h.getEntry(1) + " h z " + h.getEntry(2) + " z: r " + z.getEntry(0) + " phi " + z.getEntry(1) + " z " + z.getEntry(2) ); - + measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 + measurementMatrix = H(stateEstimation, indicator);//6x1 + h = h(stateEstimation, indicator);//1x1 + z = indicator.hit.get_Vector();//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); @@ -129,7 +116,6 @@ public void correct(Indicator indicator) { RealMatrix tmpMatrix = identity.subtract(kalmanGain.multiply(measurementMatrix)); errorCovariance = tmpMatrix.multiply(errorCovariance.multiply(tmpMatrix.transpose())).add(kalmanGain.multiply(measurementNoise.multiply(kalmanGain.transpose()))); - //System.out.println(" state after: (" + stateEstimation.getEntry(0) + ", " + stateEstimation.getEntry(1) + ", " + stateEstimation.getEntry(2) + ", " + stateEstimation.getEntry(3) + ", " + stateEstimation.getEntry(4) + ", " + stateEstimation.getEntry(5) + ");" ); // Give back to the stepper the new stateEstimation stepper.y = stateEstimation.toArray(); } @@ -139,19 +125,13 @@ public double residual(Indicator indicator) { return indicator.hit.doca()-d; } - public double wire_sign(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 + //function for left-right disambiguation + public double wire_sign(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 double phi_state = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)); double phi_wire = indicator.hit.phi(stateEstimation.getEntry(2)); - //System.out.println(" phi state " + phi_state + " phi wire " + phi_wire);// + " phi state alt? " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))); return (phi_state-phi_wire)/Math.abs(phi_state-phi_wire) ; } - // public RealMatrix wire_sign_mat(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 - // double phi_state = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)); - // double phi_wire = indicator.hit.phi(stateEstimation.getEntry(2)); - // System.out.println(" phi state " + phi_state + " phi wire " + phi_wire);// + " phi state alt? " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))); - // return MatrixUtils.createRealMatrix(new double[][]{{(phi_state-phi_wire)/Math.abs(phi_state-phi_wire)}}); - // } public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } @@ -193,63 +173,14 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { return new double[]{dxdi, dydi, dzdi, dpxdi, dpydi, dpzdi}; } - //measurement matrix in cylindrical coordinates: r, phi, z - private RealVector h(RealVector x, Indicator indicator) { - //As per my understanding: d -> r wire; phi -> phi wire, z unconstrained - double xx = x.getEntry(0); - double yy = x.getEntry(1); - return MatrixUtils.createRealVector(new double[]{Math.hypot(xx, yy), Math.atan2(yy, xx), x.getEntry(2)}); - //return MatrixUtils.createRealVector(new double[]{Math.hypot(xx, yy), Math.atan2(yy, xx)}); - } - - //measurement matrix in 1 dimension: minimize distance - doca - private RealVector h_simple(RealVector x, Indicator indicator) { + //measurement matrix in 1 dimension: minimize distance - doca + private RealVector h(RealVector x, Indicator indicator) { double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); return MatrixUtils.createRealVector(new double[]{d});//would need to have this 3x3 } - // private RealVector h_simple_sign(RealVector x, Indicator indicator, double sign) { - // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - // return MatrixUtils.createRealVector(new double[]{d*sign});//would need to have this 3x3 - // } - - //measurement noise matrix in cylindrical coordinates: r, phi, z - private RealMatrix H(RealVector x, Indicator indicator) { - // dphi/dx - 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} - }); - } - - //measurement matrix in 1 dimension: minimize distance - doca - private RealMatrix H_simple(RealVector x, Indicator indicator) { + //measurement matrix in 1 dimension: minimize distance - doca + private RealMatrix H(RealVector x, Indicator indicator) { double ddocadx = subfunctionH(x, indicator, 0); double ddocady = subfunctionH(x, indicator, 1); @@ -271,8 +202,8 @@ private RealMatrix H_simple(RealVector x, Indicator indicator) { x_plus.setEntry(i, x_plus.getEntry(i) + h); x_minus.setEntry(i, x_minus.getEntry(i) - h); - double doca_plus = h_simple(x_plus, indicator).getEntry(0); - double doca_minus = h_simple(x_minus, indicator).getEntry(0); + double doca_plus = h(x_plus, indicator).getEntry(0); + double doca_minus = h(x_minus, indicator).getEntry(0); return (doca_plus - doca_minus) / (2 * h); } 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 4663d88394..cfc759acba 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 @@ -63,16 +63,6 @@ private void propagation(ArrayList tracks, DataEvent event) { } } - - /* - Writer hitsWriter = new FileWriter("hits.dat"); - for (Point3D p : sim_hits) { - hitsWriter.write("" + p.x() + ", " + p.y() + ", " + p.z() + '\n'); - } - hitsWriter.close(); - */ - - // Initialization --------------------------------------------------------------------- final double magfield = +50; final PDGParticle proton = PDGDatabase.getParticleById(2212); @@ -94,17 +84,13 @@ private void propagation(ArrayList tracks, DataEvent event) { final double pz0 = tracks.get(0).get_pz(); final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; - //System.out.println("y = " + x0 + ", " + y0 + ", " + z0 + ", " + px0 + ", " + py0 + ", " + pz0 + "; p = " + p_init); // EPAF: *the line below is for TEST ONLY!!!* //double[] y = new double[]{vxmc, vymc, vzmc, pxmc, pymc, pzmc}; - //System.out.println("y = " + vxmc + ", " + vymc + ", " + vzmc + ", " + pxmc + ", " + pymc + ", " + pzmc + "; p = " + java.lang.Math.sqrt(pxmc*pxmc+pymc*pymc+pzmc*pzmc)); // Initialization hit - //System.out.println("tracks = " + tracks); ArrayList AHDC_hits = tracks.get(0).getHits(); ArrayList KF_hits = new ArrayList<>(); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { - //System.out.println("Superlayer = " + AHDC_hit.getSuperLayerId() + ", Layer " + AHDC_hit.getLayerId() + ", Wire " + AHDC_hit.getWireId() + ", Nwires " + AHDC_hit.getNbOfWires() + ", Radius " + AHDC_hit.getRadius() + ", DOCA " + AHDC_hit.getDoca()); Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); hit.setHitIdx(AHDC_hit.getId()); // Do delete hit with same radius @@ -117,17 +103,6 @@ private void propagation(ArrayList tracks, DataEvent event) { // if (!aleardyHaveR) KF_hits.add(hit); } - //AHDC_hits = tracks.get(0).getHits(); - - /* - Writer hitsWiresWriter = new FileWriter("hits_wires.dat"); - for (Hit h : KF_hits) { - hitsWiresWriter.write("" + h.getSuperLayer() + ", " + h.getLayer() + ", " + h.getWire() + ", " + h.getDoca() + ", " + h.getNumWires() + ", " + h.getR() + '\n'); - } - hitsWiresWriter.close(); - */ - - //System.out.println("KF_hits = " + KF_hits); final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap); @@ -145,73 +120,32 @@ private void propagation(ArrayList tracks, DataEvent event) { RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); KFitter kFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); - /* - Stepper stepper_fisrt = new Stepper(y); - Writer writer_first = new FileWriter("track_first.dat"); - for (Indicator indicator : forwardIndicators) { - stepper_fisrt.initialize(indicator); - propagator.propagateAndWrite(stepper_fisrt, indicator, writer_first); - } - writer_first.close(); - - - - System.out.println("--------- BackWard propagation !! ---------"); - - Writer writer_back = new FileWriter("track_back.dat"); - for (Indicator indicator : backwardIndicators) { - stepper.initialize(indicator); - propagator.propagateAndWrite(stepper, indicator, writer_back); - } - writer_back.close(); - */ - for (int k = 0; k < Niter; k++) { - //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: //kFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { kFitter.predict(indicator); - //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); if (indicator.haveAHit()) { - if( k==0 && indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ - //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); - AHDC_hit.setResidualPrefit(kFitter.residual(indicator)); - } + if( k==0 && indicator.hit.getHitIdx()>0){ + //sign[] = kFitter.wire_sign(indicator); + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(kFitter.residual(indicator)); + } } - //System.out.println( "AHDC hit corresponding index " + indicator.hit.getHitIdx() + " pre-fit residual " + kFitter.residual(indicator) );//AHDC_hits.get( indicator.hit.getHitIdx() ).setResidualPrefit(kFitter.residual(indicator)); - } - //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); - kFitter.correct(indicator); - //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); + kFitter.correct(indicator); } } //System.out.println("--------- BackWard propagation !! ---------"); - for (Indicator indicator : backwardIndicators) { kFitter.predict(indicator); - //System.out.println("indicator R " + indicator.R + " h " + indicator.h + "; y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); if (indicator.haveAHit()) { - //System.out.println("Superlayer = " + indicator.hit.getSuperLayer() + ", Layer " + indicator.hit.getLayer() + ", Wire " + indicator.hit.getWire() + ", Nwires " + indicator.hit.getNumWires() + ", Radius " + indicator.hit.getR() + ", DOCA " + indicator.hit.getDoca()); - kFitter.correct(indicator); - //System.out.println("y = " + kFitter.getStateEstimationVector() + " p = " + kFitter.getMomentum()); + kFitter.correct(indicator); } } } - /* - Writer writer_last = new FileWriter("track_last.dat"); - for (Indicator indicator : forwardIndicators) { - stepper.initialize(indicator); - propagator.propagateAndWrite(stepper, indicator, writer_last); - } - writer_last.close(); - */ - RealVector x_out = kFitter.getStateEstimationVector(); tracks.get(0).setPositionAndMomentumForKF(x_out); @@ -220,17 +154,12 @@ private void propagation(ArrayList tracks, DataEvent event) { kFitter.predict(indicator); if (indicator.haveAHit()) { if( indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx()){ - //System.out.println(" AHDC hit superlayer " + AHDC_hit.getSuperLayerId() ); - AHDC_hit.setResidual(kFitter.residual(indicator)); + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(kFitter.residual(indicator)); } - } } } } - - //System.out.println("y_final = " + x_out + " p_final = " + kFitter.getMomentum()); } catch (Exception e) { // e.printStackTrace(); } @@ -304,7 +233,7 @@ ArrayList backwardIndicators(ArrayList hitArrayList, HashMap= 0; i--) { backwardIndicators.add(new Indicator(hitArrayList.get(i).r(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); } - backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); + backwardIndicators.add(new Indicator(3.060, 0.1, null, false, materialHashMap.get("BONuS12Gas"))); backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); Hit hit = new Hit_beam(0, 0, 0, 0, 0, 0, 0, 0); backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); From 7aec0786536e7f77b8902369528325cfabd78e3c Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Wed, 12 Feb 2025 00:27:38 -0500 Subject: [PATCH 06/28] 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 --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 53 ++++++++++++++++++- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 16 +++--- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- 3 files changed, 63 insertions(+), 8 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index f9460b7080..73dbb12f9c 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -22,7 +22,10 @@ public class Hit implements Comparable { private final double adc; private final double numWires; private final Line3D line3D; + 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 @@ -36,6 +39,7 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do 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 @@ -106,6 +110,39 @@ 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; + + //calculate the "virtual" left and right wires accounting for the DOCA + double deltaphi = Math.sin(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 @@ -114,7 +151,7 @@ public RealVector get_Vector() { } public RealMatrix get_MeasurementNoise() { - return new Array2DRowRealMatrix(new double[][]{{0.01}}); + return new Array2DRowRealMatrix(new double[][]{{0.01}}); } public double doca() { @@ -137,6 +174,12 @@ public double distance(Point3D point3D) { return this.line3D.distance(point3D).length(); } + public double distance(Point3D point3D, int sign) { + 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(); + } + @Override public int compareTo(Hit o) { System.out.println("r = " + r + " other r = " + o.r); @@ -195,6 +238,14 @@ public int getHitIdx() { public void setHitIdx(int idx) { this.hitidx = idx; } + + public int getSign() { + return hitsign; + } + + public void setSign(int sign) { + this.hitsign = sign; + } } 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 9039911ece..927083d3c0 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 @@ -81,15 +81,14 @@ public void correct(Indicator indicator) { measurementNoise = new Array2DRowRealMatrix( new double[][]{ - {9.00, 0.0000, 0.0000}, + {0.09, 0.0000, 0.0000}, {0.00, 1e10, 0.0000}, - {0.00, 0.0000, 1e10} + {0.00, 0.0000, 1.e10} });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { - // double wire_sg = wire_sign(indicator); measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 h = h(stateEstimation, indicator);//1x1 @@ -126,10 +125,14 @@ public double residual(Indicator indicator) { } //function for left-right disambiguation - public double wire_sign(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 + public int wire_sign(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 double phi_state = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)); double phi_wire = indicator.hit.phi(stateEstimation.getEntry(2)); - return (phi_state-phi_wire)/Math.abs(phi_state-phi_wire) ; + if( (phi_state-phi_wire)/Math.abs(phi_state-phi_wire)>0 ){ + return +1; + }else{ + return -1; + } } public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ @@ -176,7 +179,8 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { //measurement matrix in 1 dimension: minimize distance - doca private RealVector h(RealVector x, Indicator indicator) { double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - return MatrixUtils.createRealVector(new double[]{d});//would need to have this 3x3 + // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); + return MatrixUtils.createRealVector(new double[]{d}); } //measurement matrix in 1 dimension: minimize distance - doca 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 cfc759acba..4035c8ddc7 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 @@ -128,7 +128,7 @@ private void propagation(ArrayList tracks, DataEvent event) { kFitter.predict(indicator); if (indicator.haveAHit()) { if( k==0 && indicator.hit.getHitIdx()>0){ - //sign[] = kFitter.wire_sign(indicator); + indicator.hit.setSign( kFitter.wire_sign(indicator) ); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(kFitter.residual(indicator)); } From 4fed2cbb27ac0bd789588f7aa87b14fb9ea42648 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Wed, 12 Feb 2025 11:12:33 -0500 Subject: [PATCH 07/28] Fix of a parameter modified by mistake. --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 4035c8ddc7..65c3babfd7 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 @@ -233,7 +233,7 @@ ArrayList backwardIndicators(ArrayList hitArrayList, HashMap= 0; i--) { backwardIndicators.add(new Indicator(hitArrayList.get(i).r(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); } - backwardIndicators.add(new Indicator(3.060, 0.1, null, false, materialHashMap.get("BONuS12Gas"))); + backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); Hit hit = new Hit_beam(0, 0, 0, 0, 0, 0, 0, 0); backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); From 189cab92d2340493204784e1ff37f4fda0e2f7c2 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Wed, 12 Feb 2025 14:18:52 -0500 Subject: [PATCH 08/28] Added a second definition of BackwardIndicators in AHDC/KalmanFilter to be able to initialize a vertex. --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) 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 65c3babfd7..a59dec7806 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 @@ -105,8 +105,8 @@ private void propagation(ArrayList tracks, DataEvent event) { } final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); - final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap); - + final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, 0); + // Start propagation Stepper stepper = new Stepper(y); RungeKutta4 RK4 = new RungeKutta4(proton, numberOfVariables, B); @@ -239,4 +239,17 @@ ArrayList backwardIndicators(ArrayList hitArrayList, HashMap backwardIndicators(ArrayList hitArrayList, HashMap materialHashMap, double vz) { + ArrayList backwardIndicators = new ArrayList<>(); + //R, h, defined in mm! + for (int i = hitArrayList.size() - 2; i >= 0; i--) { + backwardIndicators.add(new Indicator(hitArrayList.get(i).r(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); + } + backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); + backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); + Hit hit = new Hit_beam(0, 0, 0, 0, 0, 0, 0, vz); + backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); + return backwardIndicators; + } } From 8bb8f90fa6e9466a5054649ca9bb9fb2b82109c1 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Wed, 12 Feb 2025 14:41:25 -0500 Subject: [PATCH 09/28] Added a flag setDefinedVertex to AHDC/KalmanFilter and KFitter to define "hit_beam" vertex. --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 8 ++++++-- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 8 ++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) 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 927083d3c0..2b62825f98 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 @@ -17,7 +17,7 @@ public class KFitter { // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; - + private boolean isvertexdefined = false; public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator) { this.stateEstimation = initialStateEstimate; @@ -78,12 +78,14 @@ public void correct(Indicator indicator) { RealMatrix measurementMatrix; RealVector h; if (indicator.R == 0.0 && !indicator.direction) { + double z_beam_res_sq = 1.e10;//in mm + if(isvertexdefined)z_beam_res_sq = 2.25;//assuming 1.5 mm resolution measurementNoise = new Array2DRowRealMatrix( new double[][]{ {0.09, 0.0000, 0.0000}, {0.00, 1e10, 0.0000}, - {0.00, 0.0000, 1.e10} + {0.00, 0.0000, z_beam_res_sq} });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 @@ -277,4 +279,6 @@ public double getMomentum() { return Math.sqrt(stateEstimation.getEntry(3) * stateEstimation.getEntry(3) + stateEstimation.getEntry(4) * stateEstimation.getEntry(4) + stateEstimation.getEntry(5) * stateEstimation.getEntry(5)); } + public void setVertexDefined(boolean isvtxdef) {isvertexdefined = isvtxdef;} + } 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 a59dec7806..bd29424006 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 @@ -34,6 +34,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event) {propagation(tracks, event);} private final int Niter = 5; + private final boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event) { @@ -104,8 +105,10 @@ private void propagation(ArrayList tracks, DataEvent event) { KF_hits.add(hit); } + double zbeam = 0; + if(IsVtxDefined)zbeam = vxmc;//test final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); - final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, 0); + final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam); // Start propagation Stepper stepper = new Stepper(y); @@ -119,7 +122,8 @@ private void propagation(ArrayList tracks, DataEvent event) { //first 3 lines in cm^2; last 3 lines in MeV^2 RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); KFitter kFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); - + kFitter.setVertexDefined(IsVtxDefined); + for (int k = 0; k < Niter; k++) { //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: From e38f1f2d36f65a48ca45329e250d28e5320f4a9c Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Thu, 13 Feb 2025 08:19:40 -0500 Subject: [PATCH 10/28] Reset Niter and ddoca step size parameters to 10, 1.e-8 respectively. --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 6 +++--- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) 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 2b62825f98..78eedb7e7e 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 @@ -79,7 +79,7 @@ public void correct(Indicator indicator) { RealVector h; if (indicator.R == 0.0 && !indicator.direction) { double z_beam_res_sq = 1.e10;//in mm - if(isvertexdefined)z_beam_res_sq = 2.25;//assuming 1.5 mm resolution + if(isvertexdefined)z_beam_res_sq = 100.0;//assuming 10. mm resolution measurementNoise = new Array2DRowRealMatrix( new double[][]{ @@ -155,7 +155,7 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { } double[] subfunctionF(Indicator indicator, Stepper stepper1, int i) throws Exception { - double h = 1e-5;//1e-8;// in mm + double h = 1e-8;// in mm Stepper stepper_plus = new Stepper(stepper1.y); Stepper stepper_minus = new Stepper(stepper1.y); @@ -201,7 +201,7 @@ private RealMatrix H(RealVector x, Indicator indicator) { } double subfunctionH(RealVector x, Indicator indicator, int i) { - double h = 1e-5;//1e-8;// in mm + double h = 1e-8;// in mm RealVector x_plus = x.copy(); RealVector x_minus = x.copy(); 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 bd29424006..617aad45da 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 @@ -33,7 +33,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event) {propagation(tracks, event);} - private final int Niter = 5; + private final int Niter = 10; private final boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event) { From 5c642efb2c440dd5eef84e686fbd633e968344ef Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sun, 16 Feb 2025 22:48:49 -0500 Subject: [PATCH 11/28] 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. --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 7 +++- .../java/org/jlab/rec/ahdc/Hit/HitReader.java | 5 +-- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 8 +++-- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 33 +++++++++++++++---- 4 files changed, 41 insertions(+), 12 deletions(-) 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 95a202b5e6..9182365273 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 @@ -9,6 +9,7 @@ public class Hit implements Comparable { private final int layerId; private final int wireId; private final double doca; + private final double adc; private double phi; private double radius; @@ -19,12 +20,14 @@ public class Hit implements Comparable { 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; @@ -135,6 +138,8 @@ public double getY() { public double getPhi() {return phi;} + public double getADC() {return adc;} + public double getResidual() { return residual; } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java index cef052afef..91e7a1daad 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java +++ b/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 _TrueAHDCHits) { this._TrueAHDCHits = _TrueAHDCHits; } -} \ No newline at end of file +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 73dbb12f9c..5698cae87e 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -19,7 +19,7 @@ public class Hit implements Comparable { 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; private final Line3D line3D_plus; @@ -151,7 +151,7 @@ public RealVector get_Vector() { } public RealMatrix get_MeasurementNoise() { - return new Array2DRowRealMatrix(new double[][]{{0.01}}); + return new Array2DRowRealMatrix(new double[][]{{0.0225}}); } public double doca() { @@ -223,6 +223,10 @@ public double getADC() { return adc; } + public void setADC(double _adc) { + this.adc = _adc; + } + public Line3D getLine3D() { return line3D; } 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 617aad45da..775bcdd1e8 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 @@ -93,15 +93,34 @@ private void propagation(ArrayList tracks, DataEvent event) { for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); + hit.setADC(AHDC_hit.getADC()); hit.setHitIdx(AHDC_hit.getId()); // Do delete hit with same radius - // boolean aleardyHaveR = false; - // for (Hit o: KF_hits){ - // if (o.r() == hit.r()){ - // aleardyHaveR = true; - // } - // } - // if (!aleardyHaveR) + boolean aleardyHaveR = false; + for (Hit o: KF_hits){ + if (o.r() == hit.r()){ + if(hit.getADC()hit.phi()){ + o.setSign(+1); + }else{ + o.setSign(-1); + } + + }else{ + if(hit.phi()>o.phi()){ + hit.setSign(+1); + }else{ + hit.setSign(-1); + } + //remove hit + KF_hits.remove(o); + } + + } + } + if (!aleardyHaveR) KF_hits.add(hit); } From 6fe49aecac9f278d3cd3c81596015ff7f4448c46 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sun, 23 Feb 2025 22:53:06 -0500 Subject: [PATCH 12/28] 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. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 9 ++++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 47 ++++++++++--------- .../java/org/jlab/rec/ahdc/Track/Track.java | 9 ++++ .../java/org/jlab/rec/service/AHDCEngine.java | 6 ++- 4 files changed, 48 insertions(+), 23 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 5698cae87e..3261f972c5 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -150,6 +150,15 @@ public RealVector get_Vector() { return new ArrayRealVector(new double[]{this.doca}); } + //hit measurement vector in 1 dimension with sign: if sign = 0, return doca, otherwise return 0 + public RealVector get_Vector(int sign) { + if(sign == 0){ + return new ArrayRealVector(new double[]{this.doca}); + }else{ + return new ArrayRealVector(new double[]{0.0}); + } + } + public RealMatrix get_MeasurementNoise() { return new Array2DRowRealMatrix(new double[][]{{0.0225}}); } 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 775bcdd1e8..077b1b7db6 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 @@ -90,7 +90,7 @@ private void propagation(ArrayList tracks, DataEvent event) { // Initialization hit ArrayList AHDC_hits = tracks.get(0).getHits(); ArrayList KF_hits = new ArrayList<>(); - + //System.out.println(" px " + px0 + " py " + py0 +" pz " + pz0 +" vz " + z0 + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size()); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); hit.setADC(AHDC_hit.getADC()); @@ -98,34 +98,37 @@ private void propagation(ArrayList tracks, DataEvent event) { // Do delete hit with same radius boolean aleardyHaveR = false; for (Hit o: KF_hits){ - if (o.r() == hit.r()){ - if(hit.getADC()hit.phi()){ - o.setSign(+1); - }else{ - o.setSign(-1); - } - - }else{ - if(hit.phi()>o.phi()){ - hit.setSign(+1); - }else{ - hit.setSign(-1); - } - //remove hit - KF_hits.remove(o); - } - - } + // if(o.phi()>hit.phi()){ + // o.setSign(+1); + // hit.setSign(-1); + // }else{ + // hit.setSign(+1); + // o.setSign(-1); + // } + // System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); + + // }else{ + // if(hit.phi()>o.phi()){ + // hit.setSign(+1); + // }else{ + // hit.setSign(-1); + // } + // //remove hit + // KF_hits.remove(o); + // } + + } } if (!aleardyHaveR) KF_hits.add(hit); } double zbeam = 0; - if(IsVtxDefined)zbeam = vxmc;//test + if(IsVtxDefined)zbeam = vzmc;//test final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java index fee40d2b11..1e9ace5a68 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java @@ -63,6 +63,15 @@ public void setPositionAndMomentum(HelixFitObject helixFitObject) { this.pz0 = helixFitObject.get_pz(); } + public void setPositionAndMomentumVec(double[] x) { + this.x0 = x[0]; + this.y0 = x[1]; + this.z0 = x[2]; + this.px0 = x[3]; + this.py0 = x[4]; + this.pz0 = x[5]; + } + public void setPositionAndMomentumForKF(RealVector x) { this.x0_kf = x.getEntry(0); this.y0_kf = x.getEntry(1); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java index 7eb6570075..a76590de6a 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java @@ -139,7 +139,7 @@ public int compare(Hit a1, Hit a2) { } for (TrackPrediction t : predictions) { - if (t.getPrediction() > 0.5) + // if (t.getPrediction() > 0.5) AHDC_Tracks.add(new Track(t.getClusters())); } } @@ -165,6 +165,10 @@ public int compare(Hit a1, Hit a2) { HelixFitJava h = new HelixFitJava(); track.setPositionAndMomentum(h.HelixFit(nbOfPoints, szPos, 1)); + // double p = 150.0;//MeV/c + // double phi = Math.atan2(szPos[0][1], szPos[0][0]); + // double x_0[] = {0.0, 0.0, 0.0, p*Math.sin(phi), p*Math.cos(phi), 0.0}; + // track.setPositionAndMomentumVec(x_0); } // VI) Kalman Filter From 4468d66487feeabb600adc97dd69e1b67d435598 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Thu, 6 Mar 2025 06:25:54 -0500 Subject: [PATCH 13/28] Substituted call of default hit vector and hit measurement functions with hit vector and measurement functions that handle hit left/right disambiguation. --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 8 ++++---- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) 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 78eedb7e7e..300797e758 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 @@ -79,7 +79,7 @@ public void correct(Indicator indicator) { RealVector h; if (indicator.R == 0.0 && !indicator.direction) { double z_beam_res_sq = 1.e10;//in mm - if(isvertexdefined)z_beam_res_sq = 100.0;//assuming 10. mm resolution + if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 4. mm resolution measurementNoise = new Array2DRowRealMatrix( new double[][]{ @@ -94,7 +94,7 @@ public void correct(Indicator indicator) { measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 h = h(stateEstimation, indicator);//1x1 - z = indicator.hit.get_Vector();//1x1 + z = indicator.hit.get_Vector(indicator.hit.getSign());//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); @@ -180,8 +180,8 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { //measurement matrix in 1 dimension: minimize distance - doca private RealVector h(RealVector x, Indicator indicator) { - double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); + //double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); return MatrixUtils.createRealVector(new double[]{d}); } 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 077b1b7db6..4fdfa9b641 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 @@ -95,6 +95,7 @@ private void propagation(ArrayList tracks, DataEvent event) { Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); hit.setADC(AHDC_hit.getADC()); hit.setHitIdx(AHDC_hit.getId()); + hit.setSign(0); // Do delete hit with same radius boolean aleardyHaveR = false; for (Hit o: KF_hits){ @@ -154,7 +155,6 @@ private void propagation(ArrayList tracks, DataEvent event) { kFitter.predict(indicator); if (indicator.haveAHit()) { if( k==0 && indicator.hit.getHitIdx()>0){ - indicator.hit.setSign( kFitter.wire_sign(indicator) ); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(kFitter.residual(indicator)); } From 6910f04c7aa0c3ffcb015e5c0a59d60d137eb44f Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sat, 8 Mar 2025 04:49:41 -0500 Subject: [PATCH 14/28] Started to reintroduce the hit sign. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 1 + .../rec/ahdc/KalmanFilter/KalmanFilter.java | 18 ++++++++---------- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 3261f972c5..a523bbb62b 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -184,6 +184,7 @@ public double distance(Point3D point3D) { } public double distance(Point3D point3D, int sign) { + //if(sign!=0)System.out.println( " r = " + this.r + " sign ? " + sign + " distance ? sign 0: " + this.line3D.distance(point3D).length() + "; sign +: " + this.line3D_plus.distance(point3D).length() + "; sign -: " + this.line3D_minus.distance(point3D).length() ); 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(); 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 4fdfa9b641..cd98358254 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 @@ -103,15 +103,14 @@ private void propagation(ArrayList tracks, DataEvent event) { // if(hit.getADC()hit.phi()){ - // o.setSign(+1); - // hit.setSign(-1); - // }else{ - // hit.setSign(+1); - // o.setSign(-1); - // } - // System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); - + if(o.phi()>hit.phi()){ + o.setSign(-1); + hit.setSign(+1); + }else{ + hit.setSign(-1); + o.setSign(+1); + } + //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); // }else{ // if(hit.phi()>o.phi()){ // hit.setSign(+1); @@ -121,7 +120,6 @@ private void propagation(ArrayList tracks, DataEvent event) { // //remove hit // KF_hits.remove(o); // } - } } if (!aleardyHaveR) From a87bc83acc8e0b39e883c695899bbfacae5f95e0 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sun, 9 Mar 2025 20:19:21 -0400 Subject: [PATCH 15/28] Save state: back to status quo before revising sign. --- .../java/org/jlab/rec/ahdc/KalmanFilter/Hit.java | 4 +++- .../org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 3 ++- .../jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 14 ++++++++++++-- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 3261f972c5..b5204781e7 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -112,7 +112,7 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do this.line3D = wireLine; //calculate the "virtual" left and right wires accounting for the DOCA - double deltaphi = Math.sin(this.doca/R_layer); + 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 @@ -184,6 +184,8 @@ public double distance(Point3D point3D) { } public double distance(Point3D point3D, int sign) { + 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(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(); 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 300797e758..8392f32e66 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 @@ -79,7 +79,7 @@ public void correct(Indicator indicator) { RealVector h; if (indicator.R == 0.0 && !indicator.direction) { double z_beam_res_sq = 1.e10;//in mm - if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 4. mm resolution + if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 2. mm resolution measurementNoise = new Array2DRowRealMatrix( new double[][]{ @@ -93,6 +93,7 @@ public void correct(Indicator indicator) { } else { measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 + //System.out.println("h(stateEstimation):"); h = h(stateEstimation, indicator);//1x1 z = indicator.hit.get_Vector(indicator.hit.getSign());//1x1 } 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 4fdfa9b641..2749541f9f 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 @@ -83,7 +83,7 @@ private void propagation(ArrayList tracks, DataEvent event) { //final double py0 = tracks.get(0).get_py(); final double pz0 = tracks.get(0).get_pz(); - final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); + //final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // EPAF: *the line below is for TEST ONLY!!!* //double[] y = new double[]{vxmc, vymc, vzmc, pxmc, pymc, pzmc}; @@ -104,13 +104,23 @@ private void propagation(ArrayList tracks, DataEvent event) { aleardyHaveR = true; // if(o.phi()>hit.phi()){ + // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ // o.setSign(+1); // hit.setSign(-1); + // }else{ + // hit.setSign(+1); + // o.setSign(-1); + // } // }else{ + // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ // hit.setSign(+1); // o.setSign(-1); + // }else{ + // o.setSign(+1); + // hit.setSign(-1); + // } // } - // System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); + //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); // }else{ // if(hit.phi()>o.phi()){ From fd2aa169965b6c1977ccdc9366384c4560b5a839 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sun, 9 Mar 2025 23:16:43 -0400 Subject: [PATCH 16/28] Added variable measurement error for hits with sign defined, with tracks on the wrong side. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 8 +++++ .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 14 +++++--- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 36 +++++++++---------- 3 files changed, 36 insertions(+), 22 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index b5204781e7..aed03a3cdf 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -163,6 +163,14 @@ public RealMatrix get_MeasurementNoise() { return new Array2DRowRealMatrix(new double[][]{{0.0225}}); } + 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() { return doca; } 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 8392f32e66..166651fac3 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 @@ -91,11 +91,17 @@ public void correct(Indicator indicator) { h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { - measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 + boolean goodsign = true; + if(indicator.hit.getSign()!=0){ + double dphi = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))-indicator.hit.phi(stateEstimation.getEntry(2));// + if(dphi*indicator.hit.getSign()<0)goodsign = false; + } + measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 //System.out.println("h(stateEstimation):"); h = h(stateEstimation, indicator);//1x1 - z = indicator.hit.get_Vector(indicator.hit.getSign());//1x1 + //z = indicator.hit.get_Vector(indicator.hit.getSign());//1x1 + z = indicator.hit.get_Vector();//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); @@ -181,8 +187,8 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { //measurement matrix in 1 dimension: minimize distance - doca private RealVector h(RealVector x, Indicator indicator) { - //double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); + double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); return MatrixUtils.createRealVector(new double[]{d}); } 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 2749541f9f..1eb6b25b83 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 @@ -101,25 +101,25 @@ private void propagation(ArrayList tracks, DataEvent event) { for (Hit o: KF_hits){ if (o.r() == hit.r()){ // if(hit.getADC()hit.phi()){ - // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - // o.setSign(+1); - // hit.setSign(-1); - // }else{ - // hit.setSign(+1); - // o.setSign(-1); - // } - // }else{ - // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - // hit.setSign(+1); - // o.setSign(-1); - // }else{ - // o.setSign(+1); - // hit.setSign(-1); - // } - // } + if(o.phi()>hit.phi()){ + if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ + o.setSign(+1); + hit.setSign(-1); + }else{ + hit.setSign(+1); + o.setSign(-1); + } + }else{ + if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ + hit.setSign(+1); + o.setSign(-1); + }else{ + o.setSign(+1); + hit.setSign(-1); + } + } //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); // }else{ From 0857a10679da574d33a8c9234a83e8e5efdecaa2 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 10 Mar 2025 15:04:01 -0400 Subject: [PATCH 17/28] 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; --- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 4 ++- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 25 ++++++++----------- 2 files changed, 13 insertions(+), 16 deletions(-) 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 166651fac3..0b1c91fa49 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 @@ -95,8 +95,10 @@ public void correct(Indicator indicator) { if(indicator.hit.getSign()!=0){ double dphi = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))-indicator.hit.phi(stateEstimation.getEntry(2));// if(dphi*indicator.hit.getSign()<0)goodsign = false; + //System.out.println(" r " + indicator.hit.r() + " phi " + indicator.hit.phi() + " phi wire (-zl/2) " + indicator.hit.phi(-150.0) + " phi wire (0) " + indicator.hit.phi(0.0) + " phi wire (+zl/2) " + indicator.hit.phi(150.) + " phi wire (z) " + indicator.hit.phi(stateEstimation.getEntry(2)) + " phi state " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)) + " sign " + indicator.hit.getSign() + " good? " + goodsign ); } measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 + //measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 //System.out.println("h(stateEstimation):"); h = h(stateEstimation, indicator);//1x1 @@ -188,7 +190,7 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { //measurement matrix in 1 dimension: minimize distance - doca private RealVector h(RealVector x, Indicator indicator) { double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - // double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); + //double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); return MatrixUtils.createRealVector(new double[]{d}); } 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 1eb6b25b83..4163d8fe0a 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 @@ -97,10 +97,10 @@ private void propagation(ArrayList tracks, DataEvent event) { hit.setHitIdx(AHDC_hit.getId()); hit.setSign(0); // Do delete hit with same radius + boolean phi_rollover = false; boolean aleardyHaveR = false; for (Hit o: KF_hits){ if (o.r() == hit.r()){ - // if(hit.getADC()hit.phi()){ @@ -108,6 +108,7 @@ private void propagation(ArrayList tracks, DataEvent event) { o.setSign(+1); hit.setSign(-1); }else{ + phi_rollover = true; hit.setSign(+1); o.setSign(-1); } @@ -116,26 +117,20 @@ private void propagation(ArrayList tracks, DataEvent event) { hit.setSign(+1); o.setSign(-1); }else{ + phi_rollover = true; o.setSign(+1); hit.setSign(-1); } } //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); - - // }else{ - // if(hit.phi()>o.phi()){ - // hit.setSign(+1); - // }else{ - // hit.setSign(-1); - // } - // //remove hit - // KF_hits.remove(o); - // } - } } - if (!aleardyHaveR) - KF_hits.add(hit); + //if(!aleardyHaveR)KF_hits.add(hit); + if (phi_rollover){ + KF_hits.add(KF_hits.size()-1, hit); + }else{ + KF_hits.add(hit); + } } double zbeam = 0; @@ -177,7 +172,7 @@ private void propagation(ArrayList tracks, DataEvent event) { for (Indicator indicator : backwardIndicators) { kFitter.predict(indicator); if (indicator.haveAHit()) { - kFitter.correct(indicator); + kFitter.correct(indicator); } } } From a0d4aaa73fa4d81716b9cb08f003a7ce8f2369af Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Tue, 11 Mar 2025 16:45:56 -0400 Subject: [PATCH 18/28] fixed once and for all the convention for hit sign: sign >0 if phi_expected state > phi_wire --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 2 +- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 13 ++++++----- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 23 ++++++++++--------- 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index aed03a3cdf..70e8bee564 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -192,7 +192,7 @@ public double distance(Point3D point3D) { } public double distance(Point3D point3D, int sign) { - if(sign!=0) + //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(sign>0)return this.line3D_plus.distance(point3D).length(); if(sign<0)return this.line3D_minus.distance(point3D).length(); 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 0b1c91fa49..f44b892bcf 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 @@ -91,12 +91,13 @@ public void correct(Indicator indicator) { h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { - boolean goodsign = true; - if(indicator.hit.getSign()!=0){ - double dphi = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))-indicator.hit.phi(stateEstimation.getEntry(2));// - if(dphi*indicator.hit.getSign()<0)goodsign = false; - //System.out.println(" r " + indicator.hit.r() + " phi " + indicator.hit.phi() + " phi wire (-zl/2) " + indicator.hit.phi(-150.0) + " phi wire (0) " + indicator.hit.phi(0.0) + " phi wire (+zl/2) " + indicator.hit.phi(150.) + " phi wire (z) " + indicator.hit.phi(stateEstimation.getEntry(2)) + " phi state " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)) + " sign " + indicator.hit.getSign() + " good? " + goodsign ); - } + //System.out.println(" hit r " + indicator.hit.r() + " hit phi " + indicator.hit.phi() + " phi wire (-zl/2) " + indicator.hit.phi(-150.0) + " phi wire (0) " + indicator.hit.phi(0.0) + " phi wire (+zl/2) " + indicator.hit.phi(150.) + " state x " + stateEstimation.getEntry(0) + " state y " + stateEstimation.getEntry(1) + " state z " + stateEstimation.getEntry(2) ); + boolean goodsign = true; + if(indicator.hit.getSign()!=0){ + double dphi = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))-indicator.hit.phi(stateEstimation.getEntry(2)); + if(dphi*indicator.hit.getSign()<0)goodsign = false; + //System.out.println(" hit r " + indicator.hit.r() + " phi wire (z) " + indicator.hit.phi(stateEstimation.getEntry(2)) + " phi state " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)) + " sign " + indicator.hit.getSign() + " good? " + goodsign ); + } measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 //measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 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 4163d8fe0a..9622e9f9e8 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 @@ -90,36 +90,37 @@ private void propagation(ArrayList tracks, DataEvent event) { // Initialization hit ArrayList AHDC_hits = tracks.get(0).getHits(); ArrayList KF_hits = new ArrayList<>(); - //System.out.println(" px " + px0 + " py " + py0 +" pz " + pz0 +" vz " + z0 + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size()); + //System.out.println(" px " + y[3] + " py " + y[4] +" pz " + y[5] +" vz " + y[2] + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size()); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getRadius(), AHDC_hit.getDoca()); hit.setADC(AHDC_hit.getADC()); hit.setHitIdx(AHDC_hit.getId()); hit.setSign(0); + //System.out.println( " r = " + hit.r() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca() ); // Do delete hit with same radius boolean phi_rollover = false; boolean aleardyHaveR = false; for (Hit o: KF_hits){ if (o.r() == hit.r()){ aleardyHaveR = true; - + //sign+ means (phi track - phi wire) > 0 if(o.phi()>hit.phi()){ if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - o.setSign(+1); - hit.setSign(-1); + o.setSign(-1); + hit.setSign(+1); }else{ phi_rollover = true; - hit.setSign(+1); - o.setSign(-1); + hit.setSign(-1); + o.setSign(+1); } }else{ if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - hit.setSign(+1); - o.setSign(-1); + hit.setSign(-1); + o.setSign(+1); }else{ phi_rollover = true; - o.setSign(+1); - hit.setSign(-1); + o.setSign(-1); + hit.setSign(+1); } } //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); @@ -127,7 +128,7 @@ private void propagation(ArrayList tracks, DataEvent event) { } //if(!aleardyHaveR)KF_hits.add(hit); if (phi_rollover){ - KF_hits.add(KF_hits.size()-1, hit); + KF_hits.add(KF_hits.size()-1, hit); }else{ KF_hits.add(hit); } From 9bf47151d3dca0bd710e2f45a056c0a3c53966a1 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Thu, 13 Mar 2025 15:27:30 -0400 Subject: [PATCH 19/28] 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. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 70e8bee564..48d530a535 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -113,11 +113,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do //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 = -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 + 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(); @@ -128,11 +128,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do 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 = -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 + 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(); @@ -167,7 +167,7 @@ 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}}); + return new Array2DRowRealMatrix(new double[][]{{this.doca*this.doca}}); } } @@ -194,8 +194,10 @@ public double distance(Point3D point3D) { public double distance(Point3D point3D, int sign) { //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(sign>0)return this.line3D_plus.distance(point3D).length(); - if(sign<0)return this.line3D_minus.distance(point3D).length(); + if(sign>0 && sign*this.hitsign<0 + )return this.line3D_plus.distance(point3D).length(); + if(sign<0 && sign*this.hitsign<0 + )return this.line3D_minus.distance(point3D).length(); return this.line3D.distance(point3D).length(); } From 87b4394df9362d8cb3bf0fde538db85c6dd3effa Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Thu, 13 Mar 2025 15:56:57 -0400 Subject: [PATCH 20/28] 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 9bf47151d3dca0bd710e2f45a056c0a3c53966a1. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 24 +++++++++---------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 48d530a535..70e8bee564 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -113,11 +113,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do //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 = -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 + 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(); @@ -128,11 +128,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do 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 = -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 + 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(); @@ -167,7 +167,7 @@ public RealMatrix get_MeasurementNoise(boolean goodsign) { if(goodsign){ return new Array2DRowRealMatrix(new double[][]{{0.0225}}); }else{ - return new Array2DRowRealMatrix(new double[][]{{this.doca*this.doca}}); + return new Array2DRowRealMatrix(new double[][]{{2*this.doca*this.doca}}); } } @@ -194,10 +194,8 @@ public double distance(Point3D point3D) { public double distance(Point3D point3D, int sign) { //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(sign>0 && sign*this.hitsign<0 - )return this.line3D_plus.distance(point3D).length(); - if(sign<0 && sign*this.hitsign<0 - )return this.line3D_minus.distance(point3D).length(); + 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(); } From 3dfa8de8fabd2e557b4a125ca9a4265795348b2b Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Fri, 14 Mar 2025 20:13:51 -0400 Subject: [PATCH 21/28] 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. --- .../java/org/jlab/rec/ahdc/KalmanFilter/Hit.java | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 70e8bee564..327b132b91 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -113,11 +113,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do //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 = -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 + 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(); @@ -128,11 +128,11 @@ public Hit(int superLayer, int layer, int wire, int numWire, double r, double do 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 = -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 + 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(); From 8e172c1fee6b7ddc9ffeda6b5ec3ae40e7b4cf48 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Fri, 14 Mar 2025 20:30:10 -0400 Subject: [PATCH 22/28] Improved the functions to calculate hit vector: returns doca if sign is 0 or if sign is good. --- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java | 4 ++-- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 327b132b91..638c1c01c7 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -151,8 +151,8 @@ public RealVector get_Vector() { } //hit measurement vector in 1 dimension with sign: if sign = 0, return doca, otherwise return 0 - public RealVector get_Vector(int sign) { - if(sign == 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}); 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 f44b892bcf..037d263785 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 @@ -103,7 +103,7 @@ public void correct(Indicator indicator) { measurementMatrix = H(stateEstimation, indicator);//6x1 //System.out.println("h(stateEstimation):"); h = h(stateEstimation, indicator);//1x1 - //z = indicator.hit.get_Vector(indicator.hit.getSign());//1x1 + //z = indicator.hit.get_Vector(indicator.hit.getSign(), goodsign);//1x1 z = indicator.hit.get_Vector();//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); From cfc8c67d8a642c161327be60781f7deb184710be Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Fri, 14 Mar 2025 21:11:53 -0400 Subject: [PATCH 23/28] Added a hit distance function with goodsign as input, and H (measurement matrix) function with goodsign as an input. --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 8 +++-- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 36 ++++++++++++++++++- 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 638c1c01c7..a07e4299df 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -191,11 +191,13 @@ public double distance(Point3D point3D) { return this.line3D.distance(point3D).length(); } - public double distance(Point3D point3D, int sign) { + 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(sign>0)return this.line3D_plus.distance(point3D).length(); - if(sign<0)return this.line3D_minus.distance(point3D).length(); + 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(); } 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 037d263785..4199abd49a 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 @@ -189,12 +189,17 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { } //measurement matrix in 1 dimension: minimize distance - doca - private RealVector h(RealVector x, Indicator indicator) { + private RealVector h(RealVector x, Indicator indicator) { double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); //double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); return MatrixUtils.createRealVector(new double[]{d}); } + private RealVector h(RealVector x, Indicator indicator, boolean goodsign) { + double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign(), goodsign); + return MatrixUtils.createRealVector(new double[]{d}); + } + //measurement matrix in 1 dimension: minimize distance - doca private RealMatrix H(RealVector x, Indicator indicator) { @@ -224,6 +229,35 @@ private RealMatrix H(RealVector x, Indicator indicator) { return (doca_plus - doca_minus) / (2 * h); } + //measurement matrix in 1 dimension: minimize distance - doca + private RealMatrix H(RealVector x, Indicator indicator, boolean goodsign) { + + double ddocadx = subfunctionH(x, indicator, 0, goodsign); + double ddocady = subfunctionH(x, indicator, 1, goodsign); + double ddocadz = subfunctionH(x, indicator, 2, goodsign); + double ddocadpx = subfunctionH(x, indicator, 3, goodsign); + double ddocadpy = subfunctionH(x, indicator, 4, goodsign); + double ddocadpz = subfunctionH(x, indicator, 5, goodsign); + + // As per my understanding: ddocadx,y,z -> = dr/dx,y,z, etc + return MatrixUtils.createRealMatrix(new double[][]{ + {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); + } + + double subfunctionH(RealVector x, Indicator indicator, int i, boolean goodsign) { + 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, indicator, goodsign).getEntry(0); + double doca_minus = h(x_minus, indicator, goodsign).getEntry(0); + + return (doca_plus - doca_minus) / (2 * h); + } + private RealVector h_beam(RealVector x) { double xx = x.getEntry(0); From 7af95c2453de0adad524e936e340cd6489838ffb Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Fri, 14 Mar 2025 21:18:22 -0400 Subject: [PATCH 24/28] Added (commented) calls of functions with sign. --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) 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 4199abd49a..1a184bf5dc 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 @@ -98,13 +98,14 @@ public void correct(Indicator indicator) { if(dphi*indicator.hit.getSign()<0)goodsign = false; //System.out.println(" hit r " + indicator.hit.r() + " phi wire (z) " + indicator.hit.phi(stateEstimation.getEntry(2)) + " phi state " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)) + " sign " + indicator.hit.getSign() + " good? " + goodsign ); } - measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 //measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 + measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 measurementMatrix = H(stateEstimation, indicator);//6x1 - //System.out.println("h(stateEstimation):"); + //measurementMatrix = H(stateEstimation, indicator, goodsign);//6x1 h = h(stateEstimation, indicator);//1x1 - //z = indicator.hit.get_Vector(indicator.hit.getSign(), goodsign);//1x1 + //h = h(stateEstimation, indicator, goodsign);//1x1 z = indicator.hit.get_Vector();//1x1 + //z = indicator.hit.get_Vector(indicator.hit.getSign(), goodsign);//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); From 6019993d53cef8debd20e9e53547a88cd7eefbb5 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 17 Mar 2025 10:43:48 -0400 Subject: [PATCH 25/28] Added a simple handle to disable reading of MC variables. --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 50 ++++++++++--------- 1 file changed, 27 insertions(+), 23 deletions(-) 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 9622e9f9e8..68635b0ebe 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 @@ -33,35 +33,39 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event) {propagation(tracks, event);} + private final boolean IsMC = false; private final int Niter = 10; private final boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event) { try { - //If simulation read MC::Particle Bank ------------------------------------------------ - DataBank bankParticle = event.getBank("MC::Particle"); - double vxmc = bankParticle.getFloat("vx", 0)*10;//mm - double vymc = bankParticle.getFloat("vy", 0)*10;//mm - double vzmc = bankParticle.getFloat("vz", 0)*10;//mm - double pxmc = bankParticle.getFloat("px", 0)*1000;//MeV - double pymc = bankParticle.getFloat("py", 0)*1000;//MeV - double pzmc = bankParticle.getFloat("pz", 0)*1000;//MeV - double p_mc = java.lang.Math.sqrt(pxmc*pxmc+pymc*pymc+pzmc*pzmc); - //System.out.println("MC track: vz: " + vzmc*10 + " px: " + pxmc*1000 + " py: " + pymc*1000 + " pz: " + pzmc*1000 + "; p = " + p_mc*1000);//convert p to MeV, v to mm - - ArrayList sim_hits = new ArrayList<>(); - sim_hits.add(new Point3D(0, 0, vzmc)); - - DataBank bankMC = event.getBank("MC::True"); - for (int i = 0; i < bankMC.rows(); i++) { - if (bankMC.getInt("pid", i) == 2212) { - float x = bankMC.getFloat("avgX", i); - float y = bankMC.getFloat("avgY", i); - float z = bankMC.getFloat("avgZ", i); - // System.out.println("r_sim = " + Math.hypot(x, y)); - sim_hits.add(new Point3D(x, y, z)); + double vz_constraint; + if(IsMC) {//If simulation read MC::Particle Bank ------------------------------------------------ + DataBank bankParticle = event.getBank("MC::Particle"); + double vxmc = bankParticle.getFloat("vx", 0)*10;//mm + double vymc = bankParticle.getFloat("vy", 0)*10;//mm + double vzmc = bankParticle.getFloat("vz", 0)*10;//mm + double pxmc = bankParticle.getFloat("px", 0)*1000;//MeV + double pymc = bankParticle.getFloat("py", 0)*1000;//MeV + double pzmc = bankParticle.getFloat("pz", 0)*1000;//MeV + double p_mc = java.lang.Math.sqrt(pxmc*pxmc+pymc*pymc+pzmc*pzmc); + //System.out.println("MC track: vz: " + vzmc*10 + " px: " + pxmc*1000 + " py: " + pymc*1000 + " pz: " + pzmc*1000 + "; p = " + p_mc*1000);//convert p to MeV, v to mm + + ArrayList sim_hits = new ArrayList<>(); + sim_hits.add(new Point3D(0, 0, vzmc)); + + DataBank bankMC = event.getBank("MC::True"); + for (int i = 0; i < bankMC.rows(); i++) { + if (bankMC.getInt("pid", i) == 2212) { + float x = bankMC.getFloat("avgX", i); + float y = bankMC.getFloat("avgY", i); + float z = bankMC.getFloat("avgZ", i); + // System.out.println("r_sim = " + Math.hypot(x, y)); + sim_hits.add(new Point3D(x, y, z)); + } } + vz_constraint = vzmc; } // Initialization --------------------------------------------------------------------- @@ -135,7 +139,7 @@ private void propagation(ArrayList tracks, DataEvent event) { } double zbeam = 0; - if(IsVtxDefined)zbeam = vzmc;//test + if(IsVtxDefined)zbeam = vz_constraint;//test final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam); From 29cad9146da864eaee3a9065ec27d44f41cf3262 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 17 Mar 2025 14:02:08 -0400 Subject: [PATCH 26/28] Rerolled to fitting with no double hit. --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) 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 68635b0ebe..6d55b7137e 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 @@ -107,35 +107,35 @@ private void propagation(ArrayList tracks, DataEvent event) { for (Hit o: KF_hits){ if (o.r() == hit.r()){ aleardyHaveR = true; - //sign+ means (phi track - phi wire) > 0 - if(o.phi()>hit.phi()){ - if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - o.setSign(-1); - hit.setSign(+1); - }else{ - phi_rollover = true; - hit.setSign(-1); - o.setSign(+1); - } - }else{ - if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - hit.setSign(-1); - o.setSign(+1); - }else{ - phi_rollover = true; - o.setSign(-1); - hit.setSign(+1); - } - } - //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); + // //sign+ means (phi track - phi wire) > 0 + // if(o.phi()>hit.phi()){ + // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ + // o.setSign(-1); + // hit.setSign(+1); + // }else{ + // phi_rollover = true; + // hit.setSign(-1); + // o.setSign(+1); + // } + // }else{ + // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ + // hit.setSign(-1); + // o.setSign(+1); + // }else{ + // phi_rollover = true; + // o.setSign(-1); + // hit.setSign(+1); + // } + // } + // //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); } } - //if(!aleardyHaveR)KF_hits.add(hit); - if (phi_rollover){ - KF_hits.add(KF_hits.size()-1, hit); - }else{ - KF_hits.add(hit); - } + if(!aleardyHaveR)KF_hits.add(hit); + // if (phi_rollover){ + // KF_hits.add(KF_hits.size()-1, hit); + // }else{ + // KF_hits.add(hit); + // } } double zbeam = 0; From ff95394481a4b26f9fdbf090ded885041b4d9be5 Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Mon, 17 Mar 2025 14:52:29 -0400 Subject: [PATCH 27/28] Harmonized simulation flag: - one simulation flag is declared in AHDCengine and defined as false; - it is now propagated into KalmanFilter. --- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 5 ++--- .../alert/src/main/java/org/jlab/rec/service/AHDCEngine.java | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) 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 6d55b7137e..85d23fedef 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 @@ -31,13 +31,12 @@ public class KalmanFilter { - public KalmanFilter(ArrayList tracks, DataEvent event) {propagation(tracks, event);} + public KalmanFilter(ArrayList tracks, DataEvent event, boolean IsMC) {propagation(tracks, event, IsMC);} - private final boolean IsMC = false; private final int Niter = 10; private final boolean IsVtxDefined = false; - private void propagation(ArrayList tracks, DataEvent event) { + private void propagation(ArrayList tracks, DataEvent event, boolean IsMC) { try { double vz_constraint; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java index a76590de6a..ac50395dea 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java @@ -173,7 +173,7 @@ public int compare(Hit a1, Hit a2) { // VI) Kalman Filter // System.out.println("AHDC_Tracks = " + AHDC_Tracks); - KalmanFilter kalmanFitter = new KalmanFilter(AHDC_Tracks, event); + KalmanFilter kalmanFitter = new KalmanFilter(AHDC_Tracks, event, simulation); // VII) Write bank RecoBankWriter writer = new RecoBankWriter(); From 5f21a161fec9ff2feabf9f612ec2b1db3522dffd Mon Sep 17 00:00:00 2001 From: Eric Fuchey Date: Sun, 23 Mar 2025 20:56:13 -0400 Subject: [PATCH 28/28] Added a check to read MC hits in AHDC_engine. --- .../alert/src/main/java/org/jlab/rec/service/AHDCEngine.java | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java index ac50395dea..679a8445cd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/service/AHDCEngine.java @@ -79,7 +79,9 @@ public boolean processDataEvent(DataEvent event) { HitReader hitRead = new HitReader(event, simulation); ArrayList AHDC_Hits = hitRead.get_AHDCHits(); - ArrayList TrueAHDC_Hits = hitRead.get_TrueAHDCHits(); + if(simulation){ + ArrayList TrueAHDC_Hits = hitRead.get_TrueAHDCHits(); + } //System.out.println("AHDC_Hits size " + AHDC_Hits.size()); // II) Create PreCluster