From 2daa3f5b192798e526bc3848f88ce418fa72bce8 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 10 Feb 2026 15:54:25 +0100 Subject: [PATCH 01/35] add routine to clean bad hits --- .../org/jlab/service/alert/ALERTEngine.java | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index caeeba638f..f43a1cd19b 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -2,6 +2,7 @@ import java.io.File; import java.util.ArrayList; +import java.util.Iterator; import java.util.concurrent.atomic.AtomicInteger; import org.jlab.clas.reco.ReconstructionEngine; @@ -134,7 +135,8 @@ public boolean processDataEvent(DataEvent event) { if (run.get() == 0 || (run.get() != 0 && run.get() != newRun)) { run.set(newRun); } - + + /* //Do we need to read the event vx,vy,vz? //If not, this part can be moved in the initialization of the engine. double eventVx=0,eventVy=0,eventVz=0; //They should be in CM @@ -143,6 +145,7 @@ public boolean processDataEvent(DataEvent event) { float magField[] = new float[3]; swim.BfieldLab(eventVx, eventVy, eventVz, magField); this.b = Math.sqrt(Math.pow(magField[0],2) + Math.pow(magField[1],2) + Math.pow(magField[2],2)); + */ TrackProjector projector = new TrackProjector(); projector.setB(this.b); @@ -275,10 +278,35 @@ public boolean processDataEvent(DataEvent event) { PDGParticle proton = PDGDatabase.getParticleById(2212); int Niter = 40; KalmanFilter KF = new KalmanFilter(proton, Niter); + /////////////////////////////////////////////////////// // first propagation : each AHDC_tracks will be fitted /////////////////////////////////////////////////////// KF.propagation(AHDC_tracks, event, magfield, IsMC); + + ///////////////////////////////////// + /// Clean bad hits + /// ///////////////////////////////// + //System.out.println(">>>>>>>>>>>>>>>>>>>>>>>> TEST"); + double sigma = 0.5; // mm + for (Track track : AHDC_tracks) { + ArrayList AHDC_hits = track.getHits(); + Iterator it = AHDC_hits.iterator(); + while (it.hasNext()) { + Hit hit = it.next(); + //System.out.printf("> Hit : %2d %f\n", hit.getId(), hit.getResidual()); + if (Math.abs(hit.getResidual()) > 3*sigma) { + it.remove(); + } + } + } + + /////////////////////////////////////////////////////// + // second propagation : each AHDC_tracks will be fitted + /////////////////////////////////////////////////////// + KF.set_Niter(15); + KF.propagation(AHDC_tracks, event, magfield, IsMC); + ///////////////////////////////////////////// // write the AHDC::kftrack bank in the event ///////////////////////////////////////////// From 4abf9474b796792fe466a214288f236ead096450 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 10 Feb 2026 15:57:15 +0100 Subject: [PATCH 02/35] make some methods public --- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 4 ++-- 1 file changed, 2 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 272857efbd..0e77ebffb5 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 @@ -174,6 +174,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m //System.out.println("======> Kalman Filter Error"); } } - void set_Niter(int Niter) {this.Niter = Niter;} - void set_particle(PDGParticle particle) {this.particle = particle;} + public void set_Niter(int Niter) {this.Niter = Niter;} + public void set_particle(PDGParticle particle) {this.particle = particle;} } From a7bafad6e8eeedc11955324eb51ec3274bdc1ce9 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 10 Feb 2026 16:00:33 +0100 Subject: [PATCH 03/35] uncomment --- .../src/main/java/org/jlab/service/alert/ALERTEngine.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index f43a1cd19b..5c52123acf 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -136,7 +136,7 @@ public boolean processDataEvent(DataEvent event) { run.set(newRun); } - /* + //Do we need to read the event vx,vy,vz? //If not, this part can be moved in the initialization of the engine. double eventVx=0,eventVy=0,eventVz=0; //They should be in CM @@ -145,7 +145,7 @@ public boolean processDataEvent(DataEvent event) { float magField[] = new float[3]; swim.BfieldLab(eventVx, eventVy, eventVz, magField); this.b = Math.sqrt(Math.pow(magField[0],2) + Math.pow(magField[1],2) + Math.pow(magField[2],2)); - */ + TrackProjector projector = new TrackProjector(); projector.setB(this.b); From 49d13681f494e3e226596567199524c1b3545d70 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 10 Feb 2026 17:19:36 +0100 Subject: [PATCH 04/35] save error covariance matrix in Track --- .../org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 10 +++++++--- .../src/main/java/org/jlab/rec/ahdc/Track/Track.java | 12 ++++++++++++ 2 files changed, 19 insertions(+), 3 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 0e77ebffb5..7228ef6a45 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 @@ -39,10 +39,12 @@ public class KalmanFilter { private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration private double clas_alignement = -54; + private int counter = 0; // number of utilisation of the Kalman Filter public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { try { + counter++; double vz_constraint = 0; // to be linked to the electron vertex // Initialization --------------------------------------------------------------------- @@ -51,7 +53,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m final double[] B = {0.0, 0.0, magfield / 10 * tesla}; HashMap materialHashMap = MaterialMap.generateMaterials(); // Recover the vertex of the electron - if (event.hasBank("REC::Particle")) { + if (event.hasBank("REC::Particle") && !IsVtxDefined) { // as we run the KF several times, !IsVtxDefined prevent to look for the vertex again DataBank recBank = event.getBank("REC::Particle"); int row = 0; while ((!IsVtxDefined) && row < recBank.rows()) { @@ -80,7 +82,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Initialize state vector double x0 = 0.0; double y0 = 0.0; - double z0 = IsVtxDefined ? vz_constraint : track.get_Z0(); + double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); double pz0 = track.get_pz(); @@ -98,7 +100,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Initialization of the Kalman Fitter // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 RealVector initialStateEstimate = new ArrayRealVector(stepper.y); - RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); + //RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); + RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); @@ -139,6 +142,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); + track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); // Post fit propagation (no correction) to set the residuals KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); 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 e9944b7846..cceedf148f 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 @@ -1,5 +1,7 @@ package org.jlab.rec.ahdc.Track; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.rec.ahdc.AI.InterCluster; import org.jlab.rec.ahdc.Cluster.Cluster; @@ -35,6 +37,14 @@ public class Track { private double dEdx = 0; ///< deposited energy per path length (adc/mm) private double p_drift = 0; ///< momentum in the drift region (MeV) private double path = 0; ///< length of the track (mm) + // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 + RealMatrix errorCovarianceMatrix = MatrixUtils.createRealMatrix(new double[][]{ + {50 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 50 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 900 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 100 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 100 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 900}}); // AHDC::aiprediction private int predicted_ATOF_sector = -1; @@ -187,6 +197,8 @@ public void set_trackId(int _trackId) { public double get_dEdx() {return dEdx;} public double get_p_drift() {return p_drift;} public double get_path() {return path;} + public RealMatrix getErrorCovarianceMatrix() {return errorCovarianceMatrix;} + public void setErrorCovarianceMatrix(RealMatrix errorCovarianceMatrix) {this.errorCovarianceMatrix = errorCovarianceMatrix;} // AHDC::aiprediction public void set_predicted_ATOF_sector(int s) {predicted_ATOF_sector = s;} From 74cd6d297da4da828235bed0cb0edd9a848d758d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 10 Feb 2026 18:02:31 +0100 Subject: [PATCH 05/35] remove the constraint x,y = 0 when KF counter > 2 --- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 4 ++-- 1 file changed, 2 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 7228ef6a45..1067b004ac 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 @@ -80,8 +80,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Loop over tracks for (Track track : tracks) { // Initialize state vector - double x0 = 0.0; - double y0 = 0.0; + double x0 = (counter < 2) ? 0.0 : track.get_X0(); + double y0 = (counter < 2) ? 0.0 : track.get_Y0(); double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); From 39811d7bb08ea412d3d7a1d8fae83610f6ce97ab Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 13 Feb 2026 13:22:09 +0100 Subject: [PATCH 06/35] fix vz_constraint --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 33 +++++++++---------- 1 file changed, 16 insertions(+), 17 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 1067b004ac..5577d4be5f 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 @@ -16,7 +16,7 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; -//import org.apache.commons.math3.linear.RealMatrixFormat; +import org.apache.commons.math3.linear.RealMatrixFormat; /** * This is the main routine of the Kalman Filter. The fit is done by a KFitter @@ -39,13 +39,21 @@ public class KalmanFilter { private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration private double clas_alignement = -54; + double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter + // RealMatrixFormat format = + // new RealMatrixFormat( + // "[\n", "\n]", // matrix start/end + // "[", "]", // row start/end + // ",\n", // column separator + // " ; ", // row separator + // new java.text.DecimalFormat(" 0.0000;-0.0000") + // ); public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { try { counter++; - double vz_constraint = 0; // to be linked to the electron vertex // Initialization --------------------------------------------------------------------- final int numberOfVariables = 6; @@ -80,8 +88,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Loop over tracks for (Track track : tracks) { // Initialize state vector - double x0 = (counter < 2) ? 0.0 : track.get_X0(); - double y0 = (counter < 2) ? 0.0 : track.get_Y0(); + double x0 = 0.0; + double y0 = 0.0; double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); @@ -101,7 +109,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 RealVector initialStateEstimate = new ArrayRealVector(stepper.y); //RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); - RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); + // RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); + // System.out.println(">>>>>>>>>> Error matrix: start (" + counter + ")\n" + format.format(track.getErrorCovarianceMatrix())); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); @@ -128,21 +137,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } - /*RealMatrixFormat format = - new RealMatrixFormat( - "[\n", "\n]", // matrix start/end - "[", "]", // row start/end - ",\n", // column separator - " ; ", // row separator - new java.text.DecimalFormat("0.0000") - ); - System.out.println("=====> Print error matrix"); - System.out.println(format.format(TrackFitter.getErrorCovarianceMatrix()));*/ - - RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); + // System.out.println(">>>>>>>>>> Error matrix: end (" + counter + ")\n" + format.format(track.getErrorCovarianceMatrix())); + // System.out.println("> vz : " + track.get_Z0()); // Post fit propagation (no correction) to set the residuals KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); From 7b5a5d429a45a3fd068a3c6dc97086ec0cc5cad7 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 13 Feb 2026 13:31:05 +0100 Subject: [PATCH 07/35] fix error from last commit --- .../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 5577d4be5f..1ed6967b15 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 @@ -109,8 +109,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 RealVector initialStateEstimate = new ArrayRealVector(stepper.y); //RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); - // RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); // System.out.println(">>>>>>>>>> Error matrix: start (" + counter + ")\n" + format.format(track.getErrorCovarianceMatrix())); + RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); From 9c0503e4126543c74f92bb7c78c18ba1fdb01c86 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 16 Feb 2026 14:48:15 +0100 Subject: [PATCH 08/35] use atof bar hit in the kalman filter --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 5 -- .../rec/ahdc/KalmanFilter/Hit_atofBar.java | 50 ++++++++++++++++ .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 4 +- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 21 ++++++- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 14 +++++ .../org/jlab/service/alert/ALERTEngine.java | 59 ++++++++++++++----- 6 files changed, 128 insertions(+), 25 deletions(-) create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java index 008e892104..12bc3224a4 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 @@ -188,11 +188,6 @@ public RealMatrix get_MeasurementNoise() { //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } - // a signature for KalmanFilter.Hit_beam - public RealVector get_Vector_beam() { - return null; - } - public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java new file mode 100644 index 0000000000..dac76a290f --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java @@ -0,0 +1,50 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +import org.jlab.rec.ahdc.Hit.Hit; + +/** + * Same philosophy as Hit_beam + * A weird object that just wants to be considered as a Hit + * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() + * + * @author Felix Touchte Codjo + */ +public class Hit_atofBar extends Hit { + private double x,y,z; + private double r,phi; + Line3D barline; + + public Hit_atofBar(double x, double y, double z) { + super(0,0,0,0, 0, 0, -1); + this.x = x; + this.y = y; + this.z = z; + this.r = Math.sqrt(x*x + y*y); + this.phi = Math.atan2(y, x); + barline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis + } + + @Override + public RealVector get_Vector() { + return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); + } + + @Override + public Line3D getLine() { + return barline; + } + + @Override + public double distance(Point3D point3D) { + return this.barline.distance(point3D).length(); + } + + @Override + public double getRadius() { + return r; + } +} 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 fd25c2e94d..7bc3f9971f 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 @@ -7,7 +7,7 @@ import org.jlab.rec.ahdc.Hit.Hit; /** - * A weird object that just want to be considered as a Hit + * A weird object that just wants to be considered as a Hit * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() * * @author Mathieu Ouillon @@ -31,7 +31,7 @@ public Hit_beam(double x, double y , double z) { } @Override - public RealVector get_Vector_beam() { + public RealVector get_Vector() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 7263d6638f..b60c4abfe4 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 @@ -92,7 +92,8 @@ public void correct(Hit hit) { RealMatrix measurementMatrix; RealVector h; // check if the hit is the beamline - if (hit.getRadius() < 1) { + //if (hit.getRadius() < 1) { + if (hit instanceof Hit_beam) { // the diagonal elements are the squared errors in r, phi, z measurementNoise = new Array2DRowRealMatrix( @@ -103,8 +104,22 @@ public void correct(Hit hit) { });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 - z = hit.get_Vector_beam();//0! - } + z = hit.get_Vector();//0! + } + // check if the hit is an ATOF bar + else if (hit instanceof Hit_atofBar) { + // the diagonal elements are the squared errors in r, phi, z + measurementNoise = + new Array2DRowRealMatrix( + new double[][]{ + {9, 0.0000, 0.0000}, + {0.00, 9, 0.0000}, + {0.00, 0.0000, 100} + });//3x3 + measurementMatrix = H_beam(stateEstimation);//6x3 + h = h_beam(stateEstimation);//3x1 + z = hit.get_Vector();//0! + } // else, it is an AHDC hits else { measurementNoise = hit.get_MeasurementNoise();//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 1ed6967b15..92354895ef 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 @@ -15,6 +15,7 @@ import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; +import org.jlab.rec.ahdc.KalmanFilter.Hit_atofBar; import org.apache.commons.math3.linear.RealMatrixFormat; @@ -49,6 +50,7 @@ public class KalmanFilter { // " ; ", // row separator // new java.text.DecimalFormat(" 0.0000;-0.0000") // ); + HashMap ATOF_hits = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -122,6 +124,14 @@ public void propagation(ArrayList tracks, DataEvent event, final double m TrackFitter.correct(hit); } + // Forward propagation towards the ATOF bar + { + Hit hit = ATOF_hits.get(track.get_trackId()); + if (hit != null) { + TrackFitter.predict(hit, true); + TrackFitter.correct(hit); + } + } // Backward propagation (last layer to first layer) for (int i = AHDC_hits.size() - 2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); @@ -178,5 +188,9 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } } public void set_Niter(int Niter) {this.Niter = Niter;} + public int get_Niter() {return this.Niter;} public void set_particle(PDGParticle particle) {this.particle = particle;} + public PDGParticle get_particle() {return this.particle;} + public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; + public HashMap get_ATOF_hits() { return this.ATOF_hits;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 5c52123acf..4573f98074 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -4,6 +4,7 @@ import java.util.ArrayList; import java.util.Iterator; import java.util.concurrent.atomic.AtomicInteger; +import java.util.HashMap; import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.clas.swimtools.Swim; @@ -18,6 +19,7 @@ import org.jlab.rec.alert.banks.RecoBankWriter; import org.jlab.rec.alert.projections.TrackProjector; import org.jlab.rec.atof.hit.ATOFHit; +import org.jlab.rec.ahdc.KalmanFilter.Hit_atofBar; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; @@ -230,7 +232,7 @@ public boolean processDataEvent(DataEvent event) { /// Kalmam Filter /// /////////////////////////////////////// - // read the list of tracks/hits from the banks AHDC::track and AHDC::hits + /// Read the list of tracks/hits from the banks AHDC::track and AHDC::hits if (!event.hasBank("AHDC::track")) {return false;} DataBank trackBank = event.getBank("AHDC::track"); DataBank hitBank = event.getBank("AHDC::hits"); @@ -271,7 +273,7 @@ public boolean processDataEvent(DataEvent event) { AHDC_tracks.get(row).setPositionAndMomentumVec(vec); AHDC_tracks.get(row).set_trackId(trackid); } - // intialise the Kalman Filter + /// Intialise the Kalman Filter double magfieldfactor = runBank.getFloat("solenoid", 0); double magfield = 50*magfieldfactor; boolean IsMC = event.hasBank("MC::Particle"); @@ -279,37 +281,64 @@ public boolean processDataEvent(DataEvent event) { int Niter = 40; KalmanFilter KF = new KalmanFilter(proton, Niter); - /////////////////////////////////////////////////////// - // first propagation : each AHDC_tracks will be fitted - /////////////////////////////////////////////////////// + /// Add ATOF hit + HashMap ATOF_hits = new HashMap<>(); + for (Pair pair : matched_ATOF_hit_id) { + int trackid = pair.getKey(); + int atofid = pair.getKey(); + if (trackid > 0 && atofid > 0) { + // look fod the sector and the layer of this ATOF wedge + int sector = -1; + int layer = -1; + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getShort("id", row) == atofid) { + sector = bank_ATOFHits.getInt("sector", row); + layer = bank_ATOFHits.getInt("layer", row); + } + } + // now look for the bar + // with the same layer id or plus/minus 1 + if (sector > 0 && layer > 0) { + double x = 0, y = 0, z = 0; + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("component", row) == 10) { + if (bank_ATOFHits.getInt("layer", row) == layer || bank_ATOFHits.getInt("layer", row) == layer-1 || bank_ATOFHits.getInt("layer", row) == layer+1) { + x = bank_ATOFHits.getFloat("x", row); + y = bank_ATOFHits.getFloat("y", row); + z = bank_ATOFHits.getFloat("z", row); + } + } + } + if (Math.abs(x) > 0 && Math.abs(y) > 0) { + Hit_atofBar hit = new Hit_atofBar(x, y, z); + ATOF_hits.put(trackid, hit); + } + } + } + } + KF.set_ATOF_hits(ATOF_hits); + + /// First propagation : each AHDC_tracks will be fitted KF.propagation(AHDC_tracks, event, magfield, IsMC); - ///////////////////////////////////// /// Clean bad hits - /// ///////////////////////////////// - //System.out.println(">>>>>>>>>>>>>>>>>>>>>>>> TEST"); double sigma = 0.5; // mm for (Track track : AHDC_tracks) { ArrayList AHDC_hits = track.getHits(); Iterator it = AHDC_hits.iterator(); while (it.hasNext()) { Hit hit = it.next(); - //System.out.printf("> Hit : %2d %f\n", hit.getId(), hit.getResidual()); if (Math.abs(hit.getResidual()) > 3*sigma) { it.remove(); } } } - /////////////////////////////////////////////////////// - // second propagation : each AHDC_tracks will be fitted - /////////////////////////////////////////////////////// + /// Second propagation : each AHDC_tracks will be fitted KF.set_Niter(15); KF.propagation(AHDC_tracks, event, magfield, IsMC); - ///////////////////////////////////////////// - // write the AHDC::kftrack bank in the event - ///////////////////////////////////////////// + /// write the AHDC::kftrack bank in the event org.jlab.rec.ahdc.Banks.RecoBankWriter ahdc_writer = new org.jlab.rec.ahdc.Banks.RecoBankWriter(); DataBank recoKFTracksBank = ahdc_writer.fillAHDCKFTrackBank(event, AHDC_tracks); event.appendBank(recoKFTracksBank); From 83bb9fa579f770cd891e5bf0bdb6c668fd827979 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 16 Feb 2026 15:05:09 +0100 Subject: [PATCH 09/35] fix bug --- .../main/java/org/jlab/rec/ahdc/Track/Track.java | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) 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 cceedf148f..75a33d7397 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 @@ -67,13 +67,15 @@ public Track(ArrayList hitslist) { this.y0 = 0.0; this.z0 = 0.0; double p = 150.0;//MeV/c - //take first hit. - Hit hit = hitslist.get(0); - double phi = Math.atan2(hit.getY(), hit.getX()); - //hitslist. - this.px0 = p*Math.sin(phi); - this.py0 = p*Math.cos(phi); - this.pz0 = 0.0; + if (hitslist.size() > 0) { + //take first hit. + Hit hit = hitslist.get(0); + double phi = Math.atan2(hit.getY(), hit.getX()); + //hitslist. + this.px0 = p*Math.sin(phi); + this.py0 = p*Math.cos(phi); + this.pz0 = 0.0; + } } public void setPositionAndMomentum(HelixFitObject helixFitObject) { From b3091ebce2d82717bef57fdeb8be6e235664e116 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 18 Feb 2026 11:54:51 +0100 Subject: [PATCH 10/35] update noise for atof hit --- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 4 ++-- 1 file changed, 2 insertions(+), 2 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 b60c4abfe4..08d49651d2 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 @@ -112,9 +112,9 @@ else if (hit instanceof Hit_atofBar) { measurementNoise = new Array2DRowRealMatrix( new double[][]{ - {9, 0.0000, 0.0000}, + {0, 0.0000, 0.0000}, {0.00, 9, 0.0000}, - {0.00, 0.0000, 100} + {0.00, 0.0000, 6400} });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 From 44c75407edc9018ddc96eb9bd2291dedf09fc754 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 18 Feb 2026 12:02:41 +0100 Subject: [PATCH 11/35] propagate from the atof to the last ahdc layer --- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 4 ++++ 1 file changed, 4 insertions(+) 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 92354895ef..7ccecc0b38 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 @@ -131,6 +131,10 @@ public void propagation(ArrayList tracks, DataEvent event, final double m TrackFitter.predict(hit, true); TrackFitter.correct(hit); } + // propagation to the last ahdc layer + Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); + TrackFitter.predict(hhit, false); + TrackFitter.correct(hhit); } // Backward propagation (last layer to first layer) for (int i = AHDC_hits.size() - 2; i >= 0; i--) { From 6b027da82471a5f8d72bc567a54742c31d5f4063 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 18 Feb 2026 12:04:01 +0100 Subject: [PATCH 12/35] select the first bar and propagate to the surface --- .../java/org/jlab/service/alert/ALERTEngine.java | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 4573f98074..cf8743e965 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -300,12 +300,23 @@ public boolean processDataEvent(DataEvent event) { // with the same layer id or plus/minus 1 if (sector > 0 && layer > 0) { double x = 0, y = 0, z = 0; + boolean IsFound = false; for (int row = 0; row < bank_ATOFHits.rows(); row++) { + // we can have several bars for a given wegde + // for now, the vz resolution is not better for bar hits with lower or bigger energy + // let just take the first one we see if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("component", row) == 10) { - if (bank_ATOFHits.getInt("layer", row) == layer || bank_ATOFHits.getInt("layer", row) == layer-1 || bank_ATOFHits.getInt("layer", row) == layer+1) { + if (!IsFound && (bank_ATOFHits.getInt("layer", row) == layer || bank_ATOFHits.getInt("layer", row) == layer-1 || bank_ATOFHits.getInt("layer", row) == layer+1)) { x = bank_ATOFHits.getFloat("x", row); y = bank_ATOFHits.getFloat("y", row); z = bank_ATOFHits.getFloat("z", row); + // x, y correspond to the center of the bar + // let's move to the surface + double r = Math.sqrt(x*x+y*y); + double dr = 3; // width of the bar + x = (r-dr/2)*x/r; + y = (r-dr/2)*y/r; + IsFound = true; } } } From 33b49a8308fa2b433cd10a725838acc34a04f381 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 18 Feb 2026 15:04:38 +0100 Subject: [PATCH 13/35] use atof wedges instead of bars --- .../rec/ahdc/KalmanFilter/Hit_atofBar.java | 50 -------------- .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 12 ++++ .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 23 +------ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 25 ++++--- .../org/jlab/service/alert/ALERTEngine.java | 66 +++++++++---------- 5 files changed, 59 insertions(+), 117 deletions(-) delete mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java deleted file mode 100644 index dac76a290f..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_atofBar.java +++ /dev/null @@ -1,50 +0,0 @@ -package org.jlab.rec.ahdc.KalmanFilter; - -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.RealVector; -import org.jlab.geom.prim.Line3D; -import org.jlab.geom.prim.Point3D; -import org.jlab.rec.ahdc.Hit.Hit; - -/** - * Same philosophy as Hit_beam - * A weird object that just wants to be considered as a Hit - * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() - * - * @author Felix Touchte Codjo - */ -public class Hit_atofBar extends Hit { - private double x,y,z; - private double r,phi; - Line3D barline; - - public Hit_atofBar(double x, double y, double z) { - super(0,0,0,0, 0, 0, -1); - this.x = x; - this.y = y; - this.z = z; - this.r = Math.sqrt(x*x + y*y); - this.phi = Math.atan2(y, x); - barline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis - } - - @Override - public RealVector get_Vector() { - return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); - } - - @Override - public Line3D getLine() { - return barline; - } - - @Override - public double distance(Point3D point3D) { - return this.barline.distance(point3D).length(); - } - - @Override - public double getRadius() { - return r; - } -} 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 7bc3f9971f..4d2f8bc343 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 @@ -1,6 +1,8 @@ package org.jlab.rec.ahdc.KalmanFilter; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; @@ -19,6 +21,7 @@ public class Hit_beam extends Hit { private double x,y,z; private double r,phi; Line3D beamline; + RealMatrix measurementNoise = null; public Hit_beam(double x, double y , double z) { super(0,0,0,0, 0, 0, -1); @@ -34,6 +37,15 @@ public Hit_beam(double x, double y , double z) { public RealVector get_Vector() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } + + @Override + public RealMatrix get_MeasurementNoise() { + return measurementNoise; + } + + public void set_MeasurementNoise(RealMatrix measurementNoise) { + this.measurementNoise= measurementNoise; + } @Override public Line3D getLine() { 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 08d49651d2..2dce19bc86 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,32 +95,11 @@ public void correct(Hit hit) { //if (hit.getRadius() < 1) { if (hit instanceof Hit_beam) { // the diagonal elements are the squared errors in r, phi, z - measurementNoise = - new Array2DRowRealMatrix( - new double[][]{ - {vertex_resolutions[0], 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, vertex_resolutions[1]} - });//3x3 + measurementNoise = hit.get_MeasurementNoise(); measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 z = hit.get_Vector();//0! } - // check if the hit is an ATOF bar - else if (hit instanceof Hit_atofBar) { - // the diagonal elements are the squared errors in r, phi, z - measurementNoise = - new Array2DRowRealMatrix( - new double[][]{ - {0, 0.0000, 0.0000}, - {0.00, 9, 0.0000}, - {0.00, 0.0000, 6400} - });//3x3 - measurementMatrix = H_beam(stateEstimation);//6x3 - h = h_beam(stateEstimation);//3x1 - z = hit.get_Vector();//0! - } - // else, it is an AHDC hits else { measurementNoise = hit.get_MeasurementNoise();//1x1 measurementMatrix = H(stateEstimation, hit);//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 7ccecc0b38..af6f9e6cc6 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 @@ -4,6 +4,7 @@ import java.util.Collections; import java.util.HashMap; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; @@ -15,7 +16,6 @@ import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; -import org.jlab.rec.ahdc.KalmanFilter.Hit_atofBar; import org.apache.commons.math3.linear.RealMatrixFormat; @@ -50,7 +50,7 @@ public class KalmanFilter { // " ; ", // row separator // new java.text.DecimalFormat(" 0.0000;-0.0000") // ); - HashMap ATOF_hits = null; + HashMap ATOF_hits = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -130,11 +130,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m if (hit != null) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); + // Backward propagation to the last ahdc layer + Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); + TrackFitter.predict(hhit, false); + TrackFitter.correct(hhit); } - // propagation to the last ahdc layer - Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); - TrackFitter.predict(hhit, false); - TrackFitter.correct(hhit); } // Backward propagation (last layer to first layer) for (int i = AHDC_hits.size() - 2; i >= 0; i--) { @@ -144,7 +144,14 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } // Backward propagation (first layer to beamline) { - Hit hit = new Hit_beam(0, 0, vz_constraint); + Hit_beam hit = new Hit_beam(0, 0, vz_constraint); + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {vertex_resolutions[0], 0.0000, 0.0000}, + {0.00, 1e10, 0.0000}, + {0.00, 0.0000, vertex_resolutions[1]} + });//3x3; + hit.set_MeasurementNoise(measurementNoise); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } @@ -195,6 +202,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m public int get_Niter() {return this.Niter;} public void set_particle(PDGParticle particle) {this.particle = particle;} public PDGParticle get_particle() {return this.particle;} - public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; - public HashMap get_ATOF_hits() { return this.ATOF_hits;} + public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; + public HashMap get_ATOF_hits() { return this.ATOF_hits;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index cf8743e965..c5022d3b18 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -19,12 +19,14 @@ import org.jlab.rec.alert.banks.RecoBankWriter; import org.jlab.rec.alert.projections.TrackProjector; import org.jlab.rec.atof.hit.ATOFHit; -import org.jlab.rec.ahdc.KalmanFilter.Hit_atofBar; +import org.jlab.rec.ahdc.KalmanFilter.Hit_beam; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; import org.jlab.rec.ahdc.Track.Track; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; @@ -157,7 +159,6 @@ public boolean processDataEvent(DataEvent event) { /// --------------------------------------------------------------------------------------- /// Track matching using AI --------------------------------------------------------------- - if (event == null) return false; // TODO: is it useful? if (!event.hasBank("AHDC::track")) return false; DataBank bank_AHDCtracks = event.getBank("AHDC::track"); @@ -282,46 +283,39 @@ public boolean processDataEvent(DataEvent event) { KalmanFilter KF = new KalmanFilter(proton, Niter); /// Add ATOF hit - HashMap ATOF_hits = new HashMap<>(); + HashMap ATOF_hits = new HashMap<>(); for (Pair pair : matched_ATOF_hit_id) { int trackid = pair.getKey(); int atofid = pair.getKey(); if (trackid > 0 && atofid > 0) { - // look fod the sector and the layer of this ATOF wedge - int sector = -1; - int layer = -1; + // recover the wedge for (int row = 0; row < bank_ATOFHits.rows(); row++) { if (bank_ATOFHits.getShort("id", row) == atofid) { - sector = bank_ATOFHits.getInt("sector", row); - layer = bank_ATOFHits.getInt("layer", row); - } - } - // now look for the bar - // with the same layer id or plus/minus 1 - if (sector > 0 && layer > 0) { - double x = 0, y = 0, z = 0; - boolean IsFound = false; - for (int row = 0; row < bank_ATOFHits.rows(); row++) { - // we can have several bars for a given wegde - // for now, the vz resolution is not better for bar hits with lower or bigger energy - // let just take the first one we see - if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("component", row) == 10) { - if (!IsFound && (bank_ATOFHits.getInt("layer", row) == layer || bank_ATOFHits.getInt("layer", row) == layer-1 || bank_ATOFHits.getInt("layer", row) == layer+1)) { - x = bank_ATOFHits.getFloat("x", row); - y = bank_ATOFHits.getFloat("y", row); - z = bank_ATOFHits.getFloat("z", row); - // x, y correspond to the center of the bar - // let's move to the surface - double r = Math.sqrt(x*x+y*y); - double dr = 3; // width of the bar - x = (r-dr/2)*x/r; - y = (r-dr/2)*y/r; - IsFound = true; - } - } - } - if (Math.abs(x) > 0 && Math.abs(y) > 0) { - Hit_atofBar hit = new Hit_atofBar(x, y, z); + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + // A Hit_beam can be used to store the relevant information of the ATOF wedge hit + // To Do: unify Hit and Hit_beam: e.g they should implements the same interface + Hit_beam hit = new Hit_beam(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + // // error on z (from elastics data) + // double dz2 = 1400; // mm^2 + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.set_MeasurementNoise(measurementNoise); ATOF_hits.put(trackid, hit); } } From 4f26a3f4927c149fbc5758ddd76737901561f213 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 19 Feb 2026 12:27:07 +0100 Subject: [PATCH 14/35] add atof ahdc vz shift --- .../src/main/java/org/jlab/service/alert/ALERTEngine.java | 3 +++ 1 file changed, 3 insertions(+) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index c5022d3b18..50465fef7f 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -296,6 +296,7 @@ public boolean processDataEvent(DataEvent event) { double z = bank_ATOFHits.getFloat("z", row); // A Hit_beam can be used to store the relevant information of the ATOF wedge hit // To Do: unify Hit and Hit_beam: e.g they should implements the same interface + z -= 3.23; // Include the shift between AHDC and ATOF ~ 2.5 cm (what is the origin of this shift?) Hit_beam hit = new Hit_beam(x, y, z); // error on r double wedge_width = 20; //mm @@ -315,6 +316,8 @@ public boolean processDataEvent(DataEvent event) { {0.00, dphi2, 0.0000}, {0.00, 0.0000, dz2} });//3x3; + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + // we should also include the error to select the wrong wedge ! hit.set_MeasurementNoise(measurementNoise); ATOF_hits.put(trackid, hit); } From 253be6a44eca58b96bdb1bc14f149b7441c9414d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 19 Feb 2026 17:20:15 +0100 Subject: [PATCH 15/35] fix key value access issue --- .../java/org/jlab/service/alert/ALERTEngine.java | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 50465fef7f..b7d3b26185 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -260,7 +260,7 @@ public boolean processDataEvent(DataEvent event) { AHDC_hits.add(hit); } } - AHDC_tracks.add(new Track(AHDC_hits)); + Track track = new Track(AHDC_hits); // Initialise the position and the momentum using the information of the AHDC::track // position : mm // momentum : MeV @@ -271,8 +271,9 @@ public boolean processDataEvent(DataEvent event) { double py = trackBank.getFloat("py", row); double pz = trackBank.getFloat("pz", row); double[] vec = {x, y, z, px, py, pz}; - AHDC_tracks.get(row).setPositionAndMomentumVec(vec); - AHDC_tracks.get(row).set_trackId(trackid); + track.setPositionAndMomentumVec(vec); + track.set_trackId(trackid); + AHDC_tracks.add(track); } /// Intialise the Kalman Filter double magfieldfactor = runBank.getFloat("solenoid", 0); @@ -286,7 +287,7 @@ public boolean processDataEvent(DataEvent event) { HashMap ATOF_hits = new HashMap<>(); for (Pair pair : matched_ATOF_hit_id) { int trackid = pair.getKey(); - int atofid = pair.getKey(); + int atofid = pair.getValue(); if (trackid > 0 && atofid > 0) { // recover the wedge for (int row = 0; row < bank_ATOFHits.rows(); row++) { @@ -296,7 +297,7 @@ public boolean processDataEvent(DataEvent event) { double z = bank_ATOFHits.getFloat("z", row); // A Hit_beam can be used to store the relevant information of the ATOF wedge hit // To Do: unify Hit and Hit_beam: e.g they should implements the same interface - z -= 3.23; // Include the shift between AHDC and ATOF ~ 2.5 cm (what is the origin of this shift?) + z -= 32.3; // Include the shift between AHDC and ATOF ~ 2.5 cm (what is the origin of this shift?) Hit_beam hit = new Hit_beam(x, y, z); // error on r double wedge_width = 20; //mm @@ -320,6 +321,7 @@ public boolean processDataEvent(DataEvent event) { // we should also include the error to select the wrong wedge ! hit.set_MeasurementNoise(measurementNoise); ATOF_hits.put(trackid, hit); + //System.out.println("trackid : " + trackid + " atofid : " + atofid); } } } From f76628cb5c3e79560f57ff34f85b6300b9738eb5 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 19 Feb 2026 18:47:18 +0100 Subject: [PATCH 16/35] create an interface for the hit used int the Kalman Filter --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 8 ++++++- .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 8 +++---- .../org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 17 +++++++++++++++ .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 15 ++++++------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 21 +++++++++++-------- .../rec/ahdc/KalmanFilter/Propagator.java | 7 ++----- 6 files changed, 49 insertions(+), 27 deletions(-) create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java 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 12bc3224a4..ff92997e22 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,8 +9,9 @@ import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; +import org.jlab.rec.ahdc.KalmanFilter.KFHit; -public class Hit implements Comparable { +public class Hit implements Comparable, KFHit { private final int id; private final int superLayerId; @@ -106,10 +107,12 @@ public double getDoca() { return doca; } + @Override public Line3D getLine() { return wireLine; } + @Override public double getRadius() { return radius; } @@ -174,10 +177,12 @@ public void setTrackId(int _trackId) { this.trackId = _trackId; } + @Override public RealVector get_Vector() { return new ArrayRealVector(new double[]{this.doca}); } + @Override public RealMatrix get_MeasurementNoise() { double mean_error = 0.471; // mm (no difference between adc and time) double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm @@ -188,6 +193,7 @@ public RealMatrix get_MeasurementNoise() { //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } + @Override public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java index 4d2f8bc343..c9f7932b4e 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 @@ -6,17 +6,16 @@ import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; -import org.jlab.rec.ahdc.Hit.Hit; /** - * A weird object that just wants to be considered as a Hit - * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars * * @author Mathieu Ouillon * @author Éric Fuchey * @author Felix Touchte Codjo */ -public class Hit_beam extends Hit { +public class Hit_beam implements KFHit { private double x,y,z; private double r,phi; @@ -24,7 +23,6 @@ public class Hit_beam extends Hit { RealMatrix measurementNoise = null; public Hit_beam(double x, double y , double z) { - super(0,0,0,0, 0, 0, -1); this.x = x; this.y = y; this.z = z; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java new file mode 100644 index 0000000000..563bee18e5 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -0,0 +1,17 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; +/** + * All we need to use the Kalman filter + */ +public interface KFHit { + public RealVector get_Vector(); + public RealMatrix get_MeasurementNoise(); + //public void set_MeasurementNoise(RealMatrix measurementNoise); + public Line3D getLine(); + public double distance(Point3D point3D); + public double getRadius(); +} 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 2dce19bc86..387d000cff 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 @@ -39,7 +39,7 @@ public KFitter(final RealVector initialStateEstimate, final RealMatrix initialEr this.materialHashMap = materialHashMap; } - public void predict(Hit hit, boolean direction) throws Exception { + public void predict(KFHit hit, boolean direction) throws Exception { // Initialization stepper.initialize(direction); Stepper stepper1 = new Stepper(stepper.y); @@ -86,7 +86,7 @@ public void predict(Hit hit, boolean direction) throws Exception { errorCovariance = (transitionMatrix.multiply(errorCovariance.multiply(transitionMatrixT))).add(processNoise); } - public void correct(Hit hit) { + public void correct(KFHit hit) { RealVector z; RealMatrix measurementNoise; RealMatrix measurementMatrix; @@ -132,6 +132,7 @@ public void correct(Hit hit) { stepper.y = stateEstimation.toArray(); } + // specific to AHDC hits public double residual(Hit hit) { double d = hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); return hit.getDoca()-d; @@ -141,7 +142,7 @@ public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } - private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { + private RealMatrix F(KFHit hit, Stepper stepper1) throws Exception { double[] dfdx = subfunctionF(hit, stepper1, 0); double[] dfdy = subfunctionF(hit, stepper1, 1); @@ -154,7 +155,7 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { {dfdx[0], dfdy[0], dfdz[0], dfdpx[0], dfdpy[0], dfdpz[0]}, {dfdx[1], dfdy[1], dfdz[1], dfdpx[1], dfdpy[1], dfdpz[1]}, {dfdx[2], dfdy[2], dfdz[2], dfdpx[2], dfdpy[2], dfdpz[2]}, {dfdx[3], dfdy[3], dfdz[3], dfdpx[3], dfdpy[3], dfdpz[3]}, {dfdx[4], dfdy[4], dfdz[4], dfdpx[4], dfdpy[4], dfdpz[4]}, {dfdx[5], dfdy[5], dfdz[5], dfdpx[5], dfdpy[5], dfdpz[5]}}); } - double[] subfunctionF(Hit hit, Stepper stepper1, int i) throws Exception { + double[] subfunctionF(KFHit hit, Stepper stepper1, int i) throws Exception { double h = 1e-8;// in mm Stepper stepper_plus = new Stepper(stepper1.y); Stepper stepper_minus = new Stepper(stepper1.y); @@ -179,13 +180,13 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { } // Measurement matrix in 1x1 dimension: minimize distance - doca - private RealVector h(RealVector x, Hit hit) { + private RealVector h(RealVector x, KFHit hit) { double d = hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); return MatrixUtils.createRealVector(new double[]{d}); } // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) - private RealMatrix H(RealVector x, Hit hit) { + private RealMatrix H(RealVector x, KFHit hit) { double ddocadx = subfunctionH(x, hit, 0); double ddocady = subfunctionH(x, hit, 1); @@ -198,7 +199,7 @@ private RealMatrix H(RealVector x, Hit hit) { {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); } - double subfunctionH(RealVector x, Hit hit, int i) { + double subfunctionH(RealVector x, KFHit hit, int i) { 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 af6f9e6cc6..0a5b48194d 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 @@ -42,14 +42,14 @@ public class KalmanFilter { private double clas_alignement = -54; double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter - // RealMatrixFormat format = - // new RealMatrixFormat( - // "[\n", "\n]", // matrix start/end - // "[", "]", // row start/end - // ",\n", // column separator - // " ; ", // row separator - // new java.text.DecimalFormat(" 0.0000;-0.0000") - // ); + RealMatrixFormat format = + new RealMatrixFormat( + "[\n", "\n]", // matrix start/end + "[", "]", // row start/end + ",\n", // column separator + " ; ", // row separator + new java.text.DecimalFormat(" 0.0000;-0.0000") + ); HashMap ATOF_hits = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -126,7 +126,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } // Forward propagation towards the ATOF bar { - Hit hit = ATOF_hits.get(track.get_trackId()); + Hit_beam hit = ATOF_hits.get(track.get_trackId()); if (hit != null) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); @@ -134,6 +134,9 @@ public void propagation(ArrayList tracks, DataEvent event, final double m Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); TrackFitter.predict(hhit, false); TrackFitter.correct(hhit); + //System.out.println("trackid : " + track.get_trackId() + " hit : " + format.format(hit.get_MeasurementNoise())); + } else { + //System.out.println("************** No match for this track"); } } // Backward propagation (last layer to first layer) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java index e85021addb..b2b3818cb9 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java @@ -1,13 +1,10 @@ package org.jlab.rec.ahdc.KalmanFilter; -//import java.util.Arrays; - import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Point3D; import org.jlab.clas.tracking.kalmanfilter.Material; import java.util.HashMap; -import org.jlab.rec.ahdc.Hit.Hit; /** * Is responsible of the propagation @@ -27,7 +24,7 @@ public Propagator(RungeKutta4 rungeKutta4) { } // Propagate the stepper toward the next hit - void propagate(Stepper stepper, Hit hit, HashMap materialHashMap) { + void propagate(Stepper stepper, KFHit hit, HashMap materialHashMap) { // Do not allow more than 10000 steps (very critical cases) final int maxNbOfStep = 10000; @@ -161,7 +158,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } - public RealVector f(Stepper stepper, Hit hit, HashMap materialHashMap) { + public RealVector f(Stepper stepper, KFHit hit, HashMap materialHashMap) { propagate(stepper, hit, materialHashMap); return new ArrayRealVector(stepper.y); } From 860e03da154c85784bfe14e483b9cb200ccc0faa Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 20 Feb 2026 13:37:43 +0100 Subject: [PATCH 17/35] rename Hit_beam RadialKFHit --- .../java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 2 +- .../org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 10 +++++----- .../KalmanFilter/{Hit_beam.java => RadialKFHit.java} | 10 +++------- .../main/java/org/jlab/service/alert/ALERTEngine.java | 6 +++--- 4 files changed, 12 insertions(+), 16 deletions(-) rename reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/{Hit_beam.java => RadialKFHit.java} (84%) 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 387d000cff..790a10ad12 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 @@ -93,7 +93,7 @@ public void correct(KFHit hit) { RealVector h; // check if the hit is the beamline //if (hit.getRadius() < 1) { - if (hit instanceof Hit_beam) { + if (hit instanceof RadialKFHit) { // the diagonal elements are the squared errors in r, phi, z measurementNoise = hit.get_MeasurementNoise(); measurementMatrix = H_beam(stateEstimation);//6x3 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 0a5b48194d..a751a93e68 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 @@ -50,7 +50,7 @@ public class KalmanFilter { " ; ", // row separator new java.text.DecimalFormat(" 0.0000;-0.0000") ); - HashMap ATOF_hits = null; + HashMap ATOF_hits = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -126,7 +126,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } // Forward propagation towards the ATOF bar { - Hit_beam hit = ATOF_hits.get(track.get_trackId()); + RadialKFHit hit = ATOF_hits.get(track.get_trackId()); if (hit != null) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); @@ -147,7 +147,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } // Backward propagation (first layer to beamline) { - Hit_beam hit = new Hit_beam(0, 0, vz_constraint); + RadialKFHit hit = new RadialKFHit(0, 0, vz_constraint); RealMatrix measurementNoise = new Array2DRowRealMatrix( new double[][]{ {vertex_resolutions[0], 0.0000, 0.0000}, @@ -205,6 +205,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m public int get_Niter() {return this.Niter;} public void set_particle(PDGParticle particle) {this.particle = particle;} public PDGParticle get_particle() {return this.particle;} - public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; - public HashMap get_ATOF_hits() { return this.ATOF_hits;} + public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; + public HashMap get_ATOF_hits() { return this.ATOF_hits;} } 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/RadialKFHit.java similarity index 84% rename from reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java rename to reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java index c9f7932b4e..e3ab50dc27 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/RadialKFHit.java @@ -1,6 +1,5 @@ package org.jlab.rec.ahdc.KalmanFilter; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; @@ -15,16 +14,13 @@ * @author Éric Fuchey * @author Felix Touchte Codjo */ -public class Hit_beam implements KFHit { +public class RadialKFHit implements KFHit { - private double x,y,z; - private double r,phi; + private double r,phi,z; Line3D beamline; RealMatrix measurementNoise = null; - public Hit_beam(double x, double y , double z) { - this.x = x; - this.y = y; + public RadialKFHit(double x, double y , double z) { this.z = z; this.r = Math.hypot(x,y); this.phi = Math.atan2(y,x); diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index b7d3b26185..6db16b2268 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -19,7 +19,7 @@ import org.jlab.rec.alert.banks.RecoBankWriter; import org.jlab.rec.alert.projections.TrackProjector; import org.jlab.rec.atof.hit.ATOFHit; -import org.jlab.rec.ahdc.KalmanFilter.Hit_beam; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; @@ -284,7 +284,7 @@ public boolean processDataEvent(DataEvent event) { KalmanFilter KF = new KalmanFilter(proton, Niter); /// Add ATOF hit - HashMap ATOF_hits = new HashMap<>(); + HashMap ATOF_hits = new HashMap<>(); for (Pair pair : matched_ATOF_hit_id) { int trackid = pair.getKey(); int atofid = pair.getValue(); @@ -298,7 +298,7 @@ public boolean processDataEvent(DataEvent event) { // A Hit_beam can be used to store the relevant information of the ATOF wedge hit // To Do: unify Hit and Hit_beam: e.g they should implements the same interface z -= 32.3; // Include the shift between AHDC and ATOF ~ 2.5 cm (what is the origin of this shift?) - Hit_beam hit = new Hit_beam(x, y, z); + RadialKFHit hit = new RadialKFHit(x, y, z); // error on r double wedge_width = 20; //mm double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 From 9742dd8fc525255e20835c5a3f39c6c107ac8a29 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 20 Feb 2026 15:38:33 +0100 Subject: [PATCH 18/35] improve the consistency of the code --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 41 ++++++- .../org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 15 ++- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 105 +----------------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 4 +- .../rec/ahdc/KalmanFilter/RadialKFHit.java | 58 +++++++++- .../org/jlab/service/alert/ALERTEngine.java | 17 +-- 6 files changed, 117 insertions(+), 123 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 ff92997e22..d9f222aff4 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java @@ -2,6 +2,7 @@ import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.detector.calib.utils.DatabaseConstantProvider; @@ -178,12 +179,12 @@ public void setTrackId(int _trackId) { } @Override - public RealVector get_Vector() { + public RealVector MeasurementVector() { return new ArrayRealVector(new double[]{this.doca}); } @Override - public RealMatrix get_MeasurementNoise() { + public RealMatrix MeasurementNoiseMatrix() { double mean_error = 0.471; // mm (no difference between adc and time) double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm double error_on_time = (0.4423*time + 13.7215)/(0.846038*time + 31.9867); // mm @@ -193,6 +194,42 @@ public RealMatrix get_MeasurementNoise() { //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + double d = this.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + return MatrixUtils.createRealVector(new double[]{d}); + } + + // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double ddocadx = partialProjectionMatrix(x, 0); + double ddocady = partialProjectionMatrix(x, 1); + double ddocadz = partialProjectionMatrix(x, 2); + double ddocadpx = partialProjectionMatrix(x, 3); + double ddocadpy = partialProjectionMatrix(x, 4); + double ddocadpz = partialProjectionMatrix(x, 5); + + return MatrixUtils.createRealMatrix(new double[][]{ + {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); + } + + private double partialProjectionMatrix(RealVector x, int i) { + double h = 1e-8;// in mm + RealVector x_plus = x.copy(); + RealVector x_minus = x.copy(); + + x_plus.setEntry(i, x_plus.getEntry(i) + h); + x_minus.setEntry(i, x_minus.getEntry(i) - h); + + double doca_plus = this.ProjectionFunction(x_plus).getEntry(0); + double doca_minus = this.ProjectionFunction(x_minus).getEntry(0); + + return (doca_plus - doca_minus) / (2 * h); + } + @Override public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java index 563bee18e5..6884e6f011 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -5,13 +5,20 @@ import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; /** - * All we need to use the Kalman filter + * An interface to unify the hits used the Kalman Filter (e.g AHDC hits, ATOF hits, beamline) + * + * @author Felix Touchte Codjo */ public interface KFHit { - public RealVector get_Vector(); - public RealMatrix get_MeasurementNoise(); - //public void set_MeasurementNoise(RealMatrix measurementNoise); public Line3D getLine(); public double distance(Point3D point3D); public double getRadius(); + /** Return the measurement encoded in this KFHit */ + public RealVector MeasurementVector(); + /** Return the measurement noise matrix for this this KFHit */ + public RealMatrix MeasurementNoiseMatrix(); + /** Compute the measure for a given state vector */ + public RealVector ProjectionFunction(RealVector x); + /** Compute the Jacobian matrix of the {@link #ProjectionFunction(RealVector)} with respect of the components of the state vector */ + public RealMatrix ProjectionMatrix(RealVector x); } 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 790a10ad12..3067822679 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 @@ -87,25 +87,11 @@ public void predict(KFHit hit, boolean direction) throws Exception { } public void correct(KFHit hit) { - RealVector z; - RealMatrix measurementNoise; - RealMatrix measurementMatrix; - RealVector h; - // check if the hit is the beamline - //if (hit.getRadius() < 1) { - if (hit instanceof RadialKFHit) { - // the diagonal elements are the squared errors in r, phi, z - measurementNoise = hit.get_MeasurementNoise(); - measurementMatrix = H_beam(stateEstimation);//6x3 - h = h_beam(stateEstimation);//3x1 - z = hit.get_Vector();//0! - } - else { - measurementNoise = hit.get_MeasurementNoise();//1x1 - measurementMatrix = H(stateEstimation, hit);//6x1 - h = h(stateEstimation, hit);//1x1 - z = hit.get_Vector();//1x1 - } + RealVector z = hit.MeasurementVector(); + RealVector h = hit.ProjectionFunction(stateEstimation); + + RealMatrix measurementNoise = hit.MeasurementNoiseMatrix(); + RealMatrix measurementMatrix = hit.ProjectionMatrix(stateEstimation); RealMatrix measurementMatrixT = measurementMatrix.transpose(); // S = H * P(k) * H' + R @@ -179,87 +165,6 @@ private RealMatrix F(KFHit hit, Stepper stepper1) throws Exception { return new double[]{dxdi, dydi, dzdi, dpxdi, dpydi, dpzdi}; } - // Measurement matrix in 1x1 dimension: minimize distance - doca - private RealVector h(RealVector x, KFHit hit) { - double d = hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - return MatrixUtils.createRealVector(new double[]{d}); - } - - // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) - private RealMatrix H(RealVector x, KFHit hit) { - - double ddocadx = subfunctionH(x, hit, 0); - double ddocady = subfunctionH(x, hit, 1); - double ddocadz = subfunctionH(x, hit, 2); - double ddocadpx = subfunctionH(x, hit, 3); - double ddocadpy = subfunctionH(x, hit, 4); - double ddocadpz = subfunctionH(x, hit, 5); - - return MatrixUtils.createRealMatrix(new double[][]{ - {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); - } - - double subfunctionH(RealVector x, KFHit hit, int i) { - double h = 1e-8;// in mm - RealVector x_plus = x.copy(); - RealVector x_minus = x.copy(); - - x_plus.setEntry(i, x_plus.getEntry(i) + h); - x_minus.setEntry(i, x_minus.getEntry(i) - h); - - double doca_plus = h(x_plus, hit).getEntry(0); - double doca_minus = h(x_minus, hit).getEntry(0); - - return (doca_plus - doca_minus) / (2 * h); - } - - // Measurement matrix for the beamline (Hit_beam) in dimeansion 3x1 - private RealVector h_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - double zz = x.getEntry(2); - - double r = Math.hypot(xx, yy); - double phi = Math.atan2(yy, xx); - - return MatrixUtils.createRealVector(new double[]{r, phi, zz}); - } - - // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) - private RealMatrix H_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - - double drdx = (xx) / (Math.hypot(xx, yy)); - double drdy = (yy) / (Math.hypot(xx, yy)); - double drdz = 0.0; - double drdpx = 0.0; - double drdpy = 0.0; - double drdpz = 0.0; - - double dphidx = -(yy) / (xx * xx + yy * yy); - double dphidy = (xx) / (xx * xx + yy * yy); - double dphidz = 0.0; - double dphidpx = 0.0; - double dphidpy = 0.0; - double dphidpz = 0.0; - - double dzdx = 0.0; - double dzdy = 0.0; - double dzdz = 1.0; - double dzdpx = 0.0; - double dzdpy = 0.0; - double dzdpz = 0.0; - - return MatrixUtils.createRealMatrix( - new double[][]{ - {drdx, drdy, drdz, drdpx, drdpy, drdpz}, - {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, - {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} - }); - } public RealVector getStateEstimationVector() { return stateEstimation.copy(); 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 a751a93e68..14e5d7867a 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 @@ -154,7 +154,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m {0.00, 1e10, 0.0000}, {0.00, 0.0000, vertex_resolutions[1]} });//3x3; - hit.set_MeasurementNoise(measurementNoise); + hit.setMeasurementNoise(measurementNoise); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } @@ -186,7 +186,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m for (Hit hit : AHDC_hits) { sum_adc += hit.getADC(); sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0); + chi2 += Math.pow(hit.getResidual(),2)/hit.MeasurementNoiseMatrix().getEntry(0,0); } track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java index e3ab50dc27..d215194f9d 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -1,6 +1,7 @@ package org.jlab.rec.ahdc.KalmanFilter; import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Line3D; @@ -28,16 +29,16 @@ public RadialKFHit(double x, double y , double z) { } @Override - public RealVector get_Vector() { + public RealVector MeasurementVector() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } @Override - public RealMatrix get_MeasurementNoise() { + public RealMatrix MeasurementNoiseMatrix() { return measurementNoise; } - public void set_MeasurementNoise(RealMatrix measurementNoise) { + public void setMeasurementNoise(RealMatrix measurementNoise) { this.measurementNoise= measurementNoise; } @@ -55,4 +56,55 @@ public double distance(Point3D point3D) { public double getRadius() { return r; } + + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + double zz = x.getEntry(2); + + double r = Math.hypot(xx, yy); + double phi = Math.atan2(yy, xx); + + return MatrixUtils.createRealVector(new double[]{r, phi, zz}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + double dphidx = -(yy) / (xx * xx + yy * yy); + double dphidy = (xx) / (xx * xx + yy * yy); + double dphidz = 0.0; + double dphidpx = 0.0; + double dphidpy = 0.0; + double dphidpz = 0.0; + + double dzdx = 0.0; + double dzdy = 0.0; + double dzdz = 1.0; + double dzdpx = 0.0; + double dzdpy = 0.0; + double dzdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz}, + {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, + {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} + }); + } + } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 6db16b2268..3281e05710 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -283,7 +283,7 @@ public boolean processDataEvent(DataEvent event) { int Niter = 40; KalmanFilter KF = new KalmanFilter(proton, Niter); - /// Add ATOF hit + /// Add ATOF wedge hits HashMap ATOF_hits = new HashMap<>(); for (Pair pair : matched_ATOF_hit_id) { int trackid = pair.getKey(); @@ -295,9 +295,7 @@ public boolean processDataEvent(DataEvent event) { double x = bank_ATOFHits.getFloat("x", row); double y = bank_ATOFHits.getFloat("y", row); double z = bank_ATOFHits.getFloat("z", row); - // A Hit_beam can be used to store the relevant information of the ATOF wedge hit - // To Do: unify Hit and Hit_beam: e.g they should implements the same interface - z -= 32.3; // Include the shift between AHDC and ATOF ~ 2.5 cm (what is the origin of this shift?) + z -= 32.3; // there is a shift between AHDC and ATOF (still don't know why) ! RadialKFHit hit = new RadialKFHit(x, y, z); // error on r double wedge_width = 20; //mm @@ -308,8 +306,6 @@ public boolean processDataEvent(DataEvent event) { // error on z double wedge_length = 27.7; //mm double dz2 = Math.pow(wedge_length, 2)/12; - // // error on z (from elastics data) - // double dz2 = 1400; // mm^2 RealMatrix measurementNoise = new Array2DRowRealMatrix( new double[][]{ @@ -317,11 +313,8 @@ public boolean processDataEvent(DataEvent event) { {0.00, dphi2, 0.0000}, {0.00, 0.0000, dz2} });//3x3; - // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - // we should also include the error to select the wrong wedge ! - hit.set_MeasurementNoise(measurementNoise); + hit.setMeasurementNoise(measurementNoise); ATOF_hits.put(trackid, hit); - //System.out.println("trackid : " + trackid + " atofid : " + atofid); } } } @@ -331,7 +324,7 @@ public boolean processDataEvent(DataEvent event) { /// First propagation : each AHDC_tracks will be fitted KF.propagation(AHDC_tracks, event, magfield, IsMC); - /// Clean bad hits + /// Clean AHDC bad hits double sigma = 0.5; // mm for (Track track : AHDC_tracks) { ArrayList AHDC_hits = track.getHits(); @@ -359,7 +352,7 @@ public boolean processDataEvent(DataEvent event) { AHDC_hits.addAll(track.getHits()); } DataBank recoKFHitsBank = ahdc_writer.fillAHDCHitsBank(event, AHDC_hits); - event.appendBank(recoKFHitsBank); // remark: only hits assocuated to a track are saved + event.appendBank(recoKFHitsBank); // remark: only hits associated to a track are saved return true; From b6dd780d00d1798888a13f5043557dfcb1cde713 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 20 Feb 2026 15:52:08 +0100 Subject: [PATCH 19/35] remove comments or non used imports --- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 1 - .../rec/ahdc/KalmanFilter/KalmanFilter.java | 32 +++---------------- 2 files changed, 5 insertions(+), 28 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 3067822679..03bf2894c4 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java @@ -2,7 +2,6 @@ import java.util.HashMap; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; 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 14e5d7867a..47ab905113 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 @@ -6,10 +6,8 @@ import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; -import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; import org.jlab.io.base.DataBank; @@ -17,7 +15,6 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; -import org.apache.commons.math3.linear.RealMatrixFormat; /** * This is the main routine of the Kalman Filter. The fit is done by a KFitter @@ -32,8 +29,6 @@ public class KalmanFilter { public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;} - public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private PDGParticle particle; private int Niter = 40; // number of iterations for the Kalman Filter private boolean IsVtxDefined = false; // implemented but not used yet @@ -42,14 +37,6 @@ public class KalmanFilter { private double clas_alignement = -54; double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter - RealMatrixFormat format = - new RealMatrixFormat( - "[\n", "\n]", // matrix start/end - "[", "]", // row start/end - ",\n", // column separator - " ; ", // row separator - new java.text.DecimalFormat(" 0.0000;-0.0000") - ); HashMap ATOF_hits = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -70,18 +57,17 @@ public void propagation(ArrayList tracks, DataEvent event, final double m if (recBank.getInt("pid", row) == 11) { IsVtxDefined = true; vz_constraint = 10*recBank.getFloat("vz",row) - (IsMC ? 0 : clas_alignement); // mm - //////////////////////////////////////// - /// compute electron resolution here - /// it depends en p and theta - /// the fine tuning will be done later - /// //////////////////////////////////// + + // TO DO: compute electron resolution as function of p and theta + // the fine tuning will be done later //double px = recBank.getFloat("px",row); //double py = recBank.getFloat("py",row); //double pz = recBank.getFloat("pz",row); //double p = Math.sqrt(px*px+py*py+pz*pz); //double theta = Math.acos(pz/p); + vertex_resolutions[0] = 0.09; - vertex_resolutions[1] = 64;//4 + 1e10*theta + 1e10*p; + vertex_resolutions[1] = 64; } row++; } @@ -110,8 +96,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Initialization of the Kalman Fitter // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 RealVector initialStateEstimate = new ArrayRealVector(stepper.y); - //RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); - // System.out.println(">>>>>>>>>> Error matrix: start (" + counter + ")\n" + format.format(track.getErrorCovarianceMatrix())); RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); @@ -134,9 +118,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); TrackFitter.predict(hhit, false); TrackFitter.correct(hhit); - //System.out.println("trackid : " + track.get_trackId() + " hit : " + format.format(hit.get_MeasurementNoise())); - } else { - //System.out.println("************** No match for this track"); } } // Backward propagation (last layer to first layer) @@ -164,8 +145,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); - // System.out.println(">>>>>>>>>> Error matrix: end (" + counter + ")\n" + format.format(track.getErrorCovarianceMatrix())); - // System.out.println("> vz : " + track.get_Z0()); // Post fit propagation (no correction) to set the residuals KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); @@ -198,7 +177,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m }//end of loop on track candidates } catch (Exception e) { //e.printStackTrace(); - //System.out.println("======> Kalman Filter Error"); } } public void set_Niter(int Niter) {this.Niter = Niter;} From bc38c1888d73bedc4667337c0baab895b105cae8 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 24 Feb 2026 16:31:35 +0100 Subject: [PATCH 20/35] rename variable --- .../java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java index d215194f9d..7931cc51bd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -18,14 +18,14 @@ public class RadialKFHit implements KFHit { private double r,phi,z; - Line3D beamline; + Line3D line; RealMatrix measurementNoise = null; public RadialKFHit(double x, double y , double z) { this.z = z; this.r = Math.hypot(x,y); this.phi = Math.atan2(y,x); - beamline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis + line = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis } @Override @@ -44,12 +44,12 @@ public void setMeasurementNoise(RealMatrix measurementNoise) { @Override public Line3D getLine() { - return beamline; + return line; } @Override public double distance(Point3D point3D) { - return this.beamline.distance(point3D).length(); + return this.line.distance(point3D).length(); } @Override From cd7038dcb208604ffcbf46781023e1147bd6e2c4 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 24 Feb 2026 18:04:37 +0100 Subject: [PATCH 21/35] project AHDC track to the ATOF and predict wedges --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 176 ++++++++++++++++++ .../ahdc/KalmanFilter/RadialSurfaceKFHit.java | 84 +++++++++ .../org/jlab/service/alert/ALERTEngine.java | 52 +++++- 3 files changed, 311 insertions(+), 1 deletion(-) create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java 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 47ab905113..b4f16bf14a 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 @@ -3,6 +3,7 @@ import java.util.ArrayList; import java.util.Collections; import java.util.HashMap; +import java.util.Random; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; @@ -15,6 +16,13 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Component; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.geom.prim.Point3D; +import java.util.Random; + /** * This is the main routine of the Kalman Filter. The fit is done by a KFitter @@ -38,6 +46,9 @@ public class KalmanFilter { double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter HashMap ATOF_hits = null; + HashMap> ATOF_hits_predicted = new HashMap<>(); + //ArrayList ATOF_hits_predicted = new ArrayList<>(); // list of sector, layer, component (for now : only wedges) + AlertTOFDetector ATOFdet = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -154,6 +165,73 @@ public void propagation(ArrayList tracks, DataEvent event, final double m hit.setResidual(PostFitPropagator.residual(hit)); } } + // using the PostFit Propagator: check for new ATOF wegdes (fit based) + { + // The midpoint of each wedge is locate on the same radius (distance with respect to the beamline) + // We project the AHDC track in this radius + Point3D midpoint = ATOFdet.getSector(0).getSuperlayer(1).getLayer(0).getComponent(0).getMidpoint(); + double radius = midpoint.distance(new Point3D(0,0,midpoint.z())); + RadialSurfaceKFHit hit = new RadialSurfaceKFHit(radius); + // do the propagation + PostFitPropagator.predict(hit, true); + // get state vector + double[] vec = PostFitPropagator.getStateEstimationVector().toArray(); + double vx = vec[0]; + double vy = vec[1]; + double vz = vec[2]; + //System.out.printf("Initial position in ATOF : x (%.2f) y(%.2f) z (%.2f) r(%.2f)\n", vx, vy, vz, Math.sqrt(vx*vx + vy*vy)); + // translate vz with respect the ATOF center + // Warning: we corrently have a mislagnement between AHDC and ATOF that need to be investigated or fixed + vz += 32.7; // mm , temporary + int[] atof_coordinates = predict_wedge(ATOFdet, new Point3D(vx, vy, vz)); + int sector = atof_coordinates[0]; + int layer = atof_coordinates[1]; + int wedge = atof_coordinates[2]; + + // find adjacent layer and sector + int sector_plus = sector; + int sector_minus = sector; + int layer_plus = layer+1; + int layer_minus = layer-1; + if (layer == 0) { + sector_plus = sector; + sector_minus = Math.floorMod(sector-1, 15); + layer_plus = layer+1; + layer_minus = 3; + } + else if (layer == 3) { + sector_plus = Math.floorMod(sector+1, 15); + sector_minus = sector; + layer_plus = 0; + layer_minus = layer-1; + } + // Here are all the adjacents wedges (maximum 8) + //ATOF_hits_predicted.clear(); // clean first + ArrayList listOfWedges = new ArrayList<>(); + listOfWedges.add(new int[]{sector, layer, wedge}); + if (wedge-1 >= 0) { + listOfWedges.add(new int[]{sector, layer, wedge-1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge-1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge-1}); + } + if (wedge+1 <= 9) { + listOfWedges.add(new int[]{sector, layer, wedge+1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge+1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge+1}); + } + ATOF_hits_predicted.put(track.get_trackId(), listOfWedges); + + + + // if (listOfWedges != null) { + // int nn = 0; + // for (int[] res : listOfWedges) { + // nn++; + // System.out.printf("%d) Predicted wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", nn, res[0], res[1], res[2]); + // } + // } + + } // Fill track and hit bank // TO DO : s and p_drift have to be checked to be sure they represent what we want @@ -174,9 +252,13 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_dEdx(sum_adc/s); track.set_path(s); track.set_n_hits(AHDC_hits.size()); + + + }//end of loop on track candidates } catch (Exception e) { //e.printStackTrace(); + //System.out.println("Error in Kalman Filter"); } } public void set_Niter(int Niter) {this.Niter = Niter;} @@ -185,4 +267,98 @@ public void propagation(ArrayList tracks, DataEvent event, final double m public PDGParticle get_particle() {return this.particle;} public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; public HashMap get_ATOF_hits() { return this.ATOF_hits;} + public HashMap> get_ATOF_hits_predicted() {return this.ATOF_hits_predicted;} + + // Test + public static void main(String[] args) { + // ATOF detector + AlertTOFDetector atof = (new AlertTOFFactory()).createDetectorCLAS(new DatabaseConstantProvider()); + + int Npts = 1000; + int counter = 0; + int err = 0; + for (int i = 0; i < Npts; i++) { + Random rand = new Random(); + int true_sector = rand.nextInt(15); + int true_layer = rand.nextInt(4); + int true_wedge = rand.nextInt(10); + + Component comp = atof.getSector(true_sector).getSuperlayer(1).getLayer(true_layer).getComponent(true_wedge); + // top face + Point3D p0 = comp.getVolumePoint(0); + Point3D p1 = comp.getVolumePoint(1); + Point3D p2 = comp.getVolumePoint(2); + Point3D p3 = comp.getVolumePoint(3); + // bottom face + Point3D p4 = comp.getVolumePoint(4); + Point3D p5 = comp.getVolumePoint(5); + Point3D p6 = comp.getVolumePoint(6); + Point3D p7 = comp.getVolumePoint(7); + + // Random point int he current wedge volume + //Point3D pt = p0.lerp(p7, Math.random()); + //Point3D pt = comp.getMidpoint(); + double t0 = rand.nextDouble(1); + double t1 = rand.nextDouble(1-t0); + double t2 = rand.nextDouble(1-t0-t1); + double t3 = 1-t0-t1-t2; + //System.out.printf("t0 + t1 + t2 + t3 = %f\n", t0+t1+t2+t3); + double t4 = rand.nextDouble(1); + double t5 = rand.nextDouble(1-t4); + double t6 = rand.nextDouble(1-t4-t5); + double t7 = 1-t4-t5-t6; + //System.out.printf("t4 + t5 + t6 + t7 = %f\n", t4+t5+t6+t7); + double x_top = t0*p0.x() + t1*p1.x() + t2*p2.x() + t3*p3.x(); + double y_top = t0*p0.y() + t1*p1.y() + t2*p2.y() + t3*p3.y(); + double x_bot = t4*p4.x() + t5*p5.x() + t6*p6.x() + t7*p7.x(); + double y_bot = t4*p4.y() + t5*p5.y() + t6*p6.y() + t7*p7.y(); + Point3D pt_top = new Point3D(x_top, y_top, p0.z()); + Point3D pt_bot = new Point3D(x_bot, y_bot, p4.z()); + Point3D pt = pt_top.lerp(pt_bot, Math.random()); + //System.out.printf("distance from midpoint : %f\n", comp.getMidpoint().distance(pt)); + + // Test the algoritm + int[] res = KalmanFilter.predict_wedge(atof, pt); + + if (res[0] == true_sector && res[1] == true_layer && res[2] == true_wedge) { + counter++; + } else { + err++; + System.out.printf("%d) Initial wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, true_sector, true_layer, true_wedge); + System.out.printf("%d) Predicted wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, res[0], res[1], res[2]); + } + } + System.out.printf("Nb of testing : %d\n", Npts); + System.out.printf("Nb of success : %d (%.2f %%)\n", counter, 100.0*counter/Npts); + } + + static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { + // find the wedge + int wedge = -1; + double dz = 1e10; + for (int c = 0; c < atof.getSector(0).getSuperlayer(1).getLayer(0).getNumComponents(); c++) { + Point3D midpoint = atof.getSector(0).getSuperlayer(1).getLayer(0).getComponent(c).getMidpoint(); + if (Math.abs(midpoint.z()-pt.z()) < dz) { + dz = Math.abs(midpoint.z()-pt.z()); + wedge = c; + } + } + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(1).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(1).getLayer(l).getComponent(wedge).getMidpoint(); + if (midpoint.distance(pt) < d) { + d = midpoint.distance(pt); + sector = s; + layer = l; + } + } + } + return new int[] {sector, layer, wedge}; + } + + public void set_ATOF_detector(AlertTOFDetector atof) { this.ATOFdet = atof;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java new file mode 100644 index 0000000000..c3e026cec5 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java @@ -0,0 +1,84 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; + +/** + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ +public class RadialSurfaceKFHit implements KFHit { + + private double r; + + public RadialSurfaceKFHit(double r) { + this.r = r; + } + + @Override + public RealVector MeasurementVector() { + return new ArrayRealVector(new double[] {this.r}); + } + + @Override + public RealMatrix MeasurementNoiseMatrix() { + return new Array2DRowRealMatrix(new double[][]{{1e-8}}); + } + + @Override + public Line3D getLine() { + return null; + } + + @Override + public double distance(Point3D point3D) { + return Math.abs(this.r - Math.hypot(point3D.x(), point3D.y())); + } + + @Override + public double getRadius() { + return r; + } + + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double r = Math.hypot(xx, yy); + + return MatrixUtils.createRealVector(new double[]{r}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz} + }); + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 3281e05710..9e964466dd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -11,6 +11,7 @@ import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.base.Detector; import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.io.hipo.HipoDataSource; @@ -57,7 +58,7 @@ public class ALERTEngine extends ReconstructionEngine { */ private RecoBankWriter rbc; - Detector ATOF; // ALERT ATOF detector + private AlertTOFDetector ATOF; // ALERT ATOF detector private AlertDCDetector AHDC; // ALERT AHDC detector /** @@ -320,10 +321,59 @@ public boolean processDataEvent(DataEvent event) { } } KF.set_ATOF_hits(ATOF_hits); + KF.set_ATOF_detector(ATOF); /// First propagation : each AHDC_tracks will be fitted KF.propagation(AHDC_tracks, event, magfield, IsMC); + // Get ATOF wedges predicted + HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); + for (Track track : AHDC_tracks) { + int trackid = track.get_trackId(); + // if the AI has no match for this track, look for the predicted ones + if (ATOF_hits.get(trackid) == null) { + ArrayList predicted_wedges = ATOF_hits_predicted.get(trackid); + boolean IsHitSelected = false; + int i = 0; + while (i < predicted_wedges.size() && !IsHitSelected) { + int sector = predicted_wedges.get(i)[0]; + int layer = predicted_wedges.get(i)[1]; + int wedge = predicted_wedges.get(i)[2]; + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("layer", row) == layer && bank_ATOFHits.getInt("component", row) == wedge) { + // create a RadialKFHit + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + z -= 32.3; // there is a shift between AHDC and ATOF (still don't know why) ! + RadialKFHit hit = new RadialKFHit(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.setMeasurementNoise(measurementNoise); + ATOF_hits.put(trackid, hit); + + IsHitSelected = true; + break; + } + } // end loop over ATOF bank rows + } //end loop over predicted ATOF hits + } // end if + } // end loop over tracks + /// Clean AHDC bad hits double sigma = 0.5; // mm for (Track track : AHDC_tracks) { From b5643a20e7f4c6f5241c4273b3df6fe32293739d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 15:06:17 +0100 Subject: [PATCH 22/35] start projection on the 3 surfaces of the ATOF --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 183 +++++++++++------- .../org/jlab/service/alert/ALERTEngine.java | 92 ++++----- 2 files changed, 159 insertions(+), 116 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 b4f16bf14a..60b0af550d 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 @@ -41,8 +41,9 @@ public class KalmanFilter { private int Niter = 40; // number of iterations for the Kalman Filter private boolean IsVtxDefined = false; // implemented but not used yet private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} - // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration - private double clas_alignement = -54; + // mm, they the misalignement with respect to the AHDC + private double clas_alignement = +54; + private double atof_alignement = -32.7; double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter HashMap ATOF_hits = null; @@ -67,7 +68,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m while ((!IsVtxDefined) && row < recBank.rows()) { if (recBank.getInt("pid", row) == 11) { IsVtxDefined = true; - vz_constraint = 10*recBank.getFloat("vz",row) - (IsMC ? 0 : clas_alignement); // mm + vz_constraint = 10*recBank.getFloat("vz",row) + (IsMC ? 0 : clas_alignement); // mm // TO DO: compute electron resolution as function of p and theta // the fine tuning will be done later @@ -159,78 +160,57 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Post fit propagation (no correction) to set the residuals KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); + // Projection towards AHDC hits for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits hit.setResidual(PostFitPropagator.residual(hit)); } } - // using the PostFit Propagator: check for new ATOF wegdes (fit based) + // Projection towards the ATOF surfaces + // R1 : lower surface of an ATOF bar + // R2 : upper surface of an ATOF bar = lower surface of an ATOF wedge + // R3 : upper surface of an ATOF wedge + Point3D pR1 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(0); + Point3D pR2 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(2); + Point3D pR3 = ATOFdet.getSector(0).getSuperlayer(1).getLayer(0).getComponent(0).getVolumePoint(2); + double R1 = Math.hypot(pR1.x(), pR1.y()); + double R2 = Math.hypot(pR2.x(), pR2.y()); + double R3 = Math.hypot(pR3.x(), pR3.y()); + + System.out.println("ATOF surface R2 : " + R2); + System.out.println("ATOF surface R3 : " + R3); { - // The midpoint of each wedge is locate on the same radius (distance with respect to the beamline) - // We project the AHDC track in this radius - Point3D midpoint = ATOFdet.getSector(0).getSuperlayer(1).getLayer(0).getComponent(0).getMidpoint(); - double radius = midpoint.distance(new Point3D(0,0,midpoint.z())); - RadialSurfaceKFHit hit = new RadialSurfaceKFHit(radius); - // do the propagation - PostFitPropagator.predict(hit, true); - // get state vector - double[] vec = PostFitPropagator.getStateEstimationVector().toArray(); - double vx = vec[0]; - double vy = vec[1]; - double vz = vec[2]; - //System.out.printf("Initial position in ATOF : x (%.2f) y(%.2f) z (%.2f) r(%.2f)\n", vx, vy, vz, Math.sqrt(vx*vx + vy*vy)); - // translate vz with respect the ATOF center - // Warning: we corrently have a mislagnement between AHDC and ATOF that need to be investigated or fixed - vz += 32.7; // mm , temporary - int[] atof_coordinates = predict_wedge(ATOFdet, new Point3D(vx, vy, vz)); - int sector = atof_coordinates[0]; - int layer = atof_coordinates[1]; - int wedge = atof_coordinates[2]; - - // find adjacent layer and sector - int sector_plus = sector; - int sector_minus = sector; - int layer_plus = layer+1; - int layer_minus = layer-1; - if (layer == 0) { - sector_plus = sector; - sector_minus = Math.floorMod(sector-1, 15); - layer_plus = layer+1; - layer_minus = 3; - } - else if (layer == 3) { - sector_plus = Math.floorMod(sector+1, 15); - sector_minus = sector; - layer_plus = 0; - layer_minus = layer-1; - } - // Here are all the adjacents wedges (maximum 8) - //ATOF_hits_predicted.clear(); // clean first - ArrayList listOfWedges = new ArrayList<>(); - listOfWedges.add(new int[]{sector, layer, wedge}); - if (wedge-1 >= 0) { - listOfWedges.add(new int[]{sector, layer, wedge-1}); - listOfWedges.add(new int[]{sector_plus, layer_plus, wedge-1}); - listOfWedges.add(new int[]{sector_minus, layer_minus, wedge-1}); - } - if (wedge+1 <= 9) { - listOfWedges.add(new int[]{sector, layer, wedge+1}); - listOfWedges.add(new int[]{sector_plus, layer_plus, wedge+1}); - listOfWedges.add(new int[]{sector_minus, layer_minus, wedge+1}); - } - ATOF_hits_predicted.put(track.get_trackId(), listOfWedges); - - - - // if (listOfWedges != null) { - // int nn = 0; - // for (int[] res : listOfWedges) { - // nn++; - // System.out.printf("%d) Predicted wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", nn, res[0], res[1], res[2]); - // } - // } - + // From last AHDC hit to surface R1 + RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); + PostFitPropagator.predict(hitR1, true); + double[] vecR1 = PostFitPropagator.getStateEstimationVector().toArray(); + Point3D posR1 = new Point3D(vecR1[0], vecR1[1], vecR1[2]); + posR1.translateXYZ(0, 0, atof_alignement); + int[] idR1 = predict_bar(ATOFdet, posR1); + System.out.println("ATOF surface R1 : " + R1); + System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR1.x(), posR1.y(), posR1.z()); + System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR1[0], idR1[1], idR1[2]); + // From surface R1 to surface R2 + RadialSurfaceKFHit hitR2 = new RadialSurfaceKFHit(R2); + PostFitPropagator.predict(hitR2, true); + double[] vecR2 = PostFitPropagator.getStateEstimationVector().toArray(); + Point3D posR2 = new Point3D(vecR2[0], vecR2[1], vecR2[2]); + posR2.translateXYZ(0, 0, atof_alignement); + int[] idR2 = predict_wedge(ATOFdet, posR2); + System.out.println("ATOF surface R2 : " + R2); + System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR2.x(), posR2.y(), posR2.z()); + System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR2[0], idR2[1], idR2[2]); + // From surface R2 to surface R3 + RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); + PostFitPropagator.predict(hitR3, true); + double[] vecR3 = PostFitPropagator.getStateEstimationVector().toArray(); + Point3D posR3 = new Point3D(vecR3[0], vecR3[1], vecR3[2]); + posR3.translateXYZ(0, 0, atof_alignement); + int[] idR3 = predict_wedge(ATOFdet, posR3); + System.out.println("ATOF surface R3 : " + R3); + System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR3.x(), posR3.y(), posR3.z()); + System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR3[0], idR3[1], idR3[2]); } // Fill track and hit bank @@ -257,7 +237,7 @@ else if (layer == 3) { }//end of loop on track candidates } catch (Exception e) { - //e.printStackTrace(); + e.printStackTrace(); //System.out.println("Error in Kalman Filter"); } } @@ -331,7 +311,9 @@ public static void main(String[] args) { System.out.printf("Nb of testing : %d\n", Npts); System.out.printf("Nb of success : %d (%.2f %%)\n", counter, 100.0*counter/Npts); } - + /** + * @param pt is defined in the center of the ATOF + */ static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { // find the wedge int wedge = -1; @@ -360,5 +342,66 @@ static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { return new int[] {sector, layer, wedge}; } + /** + * @param pt is defined in the center of the ATOF + */ + static public int[] predict_bar(AlertTOFDetector atof, Point3D pt) { + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(0).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(0).getLayer(l).getComponent(10).getMidpoint(); + double distance = midpoint.vectorTo(pt).rho(); + if (distance < d) { + d = distance; + sector = s; + layer = l; + } + } + } + return new int[] {sector, layer, 10}; + } + + public ArrayList get_adjacent_wedges(int[] identifiers) { + int sector = identifiers[0]; + int layer = identifiers[1]; + int wedge = identifiers[2]; + + // find adjacent layer and sector + int sector_plus = sector; + int sector_minus = sector; + int layer_plus = layer+1; + int layer_minus = layer-1; + if (layer == 0) { + sector_plus = sector; + sector_minus = Math.floorMod(sector-1, 15); + layer_plus = layer+1; + layer_minus = 3; + } + else if (layer == 3) { + sector_plus = Math.floorMod(sector+1, 15); + sector_minus = sector; + layer_plus = 0; + layer_minus = layer-1; + } + // Here are all the adjacents wedges (maximum 8) + ArrayList listOfWedges = new ArrayList<>(); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge}); + if (wedge-1 >= 0) { + listOfWedges.add(new int[]{sector, layer, wedge-1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge-1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge-1}); + } + if (wedge+1 <= 9) { + listOfWedges.add(new int[]{sector, layer, wedge+1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge+1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge+1}); + } + return listOfWedges; + } + public void set_ATOF_detector(AlertTOFDetector atof) { this.ATOFdet = atof;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 9e964466dd..48dea9e98d 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -326,53 +326,53 @@ public boolean processDataEvent(DataEvent event) { /// First propagation : each AHDC_tracks will be fitted KF.propagation(AHDC_tracks, event, magfield, IsMC); - // Get ATOF wedges predicted - HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); - for (Track track : AHDC_tracks) { - int trackid = track.get_trackId(); - // if the AI has no match for this track, look for the predicted ones - if (ATOF_hits.get(trackid) == null) { - ArrayList predicted_wedges = ATOF_hits_predicted.get(trackid); - boolean IsHitSelected = false; - int i = 0; - while (i < predicted_wedges.size() && !IsHitSelected) { - int sector = predicted_wedges.get(i)[0]; - int layer = predicted_wedges.get(i)[1]; - int wedge = predicted_wedges.get(i)[2]; - for (int row = 0; row < bank_ATOFHits.rows(); row++) { - if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("layer", row) == layer && bank_ATOFHits.getInt("component", row) == wedge) { - // create a RadialKFHit - double x = bank_ATOFHits.getFloat("x", row); - double y = bank_ATOFHits.getFloat("y", row); - double z = bank_ATOFHits.getFloat("z", row); - z -= 32.3; // there is a shift between AHDC and ATOF (still don't know why) ! - RadialKFHit hit = new RadialKFHit(x, y, z); - // error on r - double wedge_width = 20; //mm - double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 - // error on phi - double open_angle = Math.toRadians(6); // deg - double dphi2 = Math.pow(open_angle, 2)/12; - // error on z - double wedge_length = 27.7; //mm - double dz2 = Math.pow(wedge_length, 2)/12; + // // Get ATOF wedges predicted + // HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); + // for (Track track : AHDC_tracks) { + // int trackid = track.get_trackId(); + // // if the AI has no match for this track, look for the predicted ones + // if (ATOF_hits.get(trackid) == null) { + // ArrayList predicted_wedges = ATOF_hits_predicted.get(trackid); + // boolean IsHitSelected = false; + // int i = 0; + // while (i < predicted_wedges.size() && !IsHitSelected) { + // int sector = predicted_wedges.get(i)[0]; + // int layer = predicted_wedges.get(i)[1]; + // int wedge = predicted_wedges.get(i)[2]; + // for (int row = 0; row < bank_ATOFHits.rows(); row++) { + // if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("layer", row) == layer && bank_ATOFHits.getInt("component", row) == wedge) { + // // create a RadialKFHit + // double x = bank_ATOFHits.getFloat("x", row); + // double y = bank_ATOFHits.getFloat("y", row); + // double z = bank_ATOFHits.getFloat("z", row); + // z -= 32.3; // there is a shift between AHDC and ATOF (still don't know why) ! + // RadialKFHit hit = new RadialKFHit(x, y, z); + // // error on r + // double wedge_width = 20; //mm + // double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // // error on phi + // double open_angle = Math.toRadians(6); // deg + // double dphi2 = Math.pow(open_angle, 2)/12; + // // error on z + // double wedge_length = 27.7; //mm + // double dz2 = Math.pow(wedge_length, 2)/12; - RealMatrix measurementNoise = new Array2DRowRealMatrix( - new double[][]{ - {dr2, 0.0000, 0.0000}, - {0.00, dphi2, 0.0000}, - {0.00, 0.0000, dz2} - });//3x3; - hit.setMeasurementNoise(measurementNoise); - ATOF_hits.put(trackid, hit); - - IsHitSelected = true; - break; - } - } // end loop over ATOF bank rows - } //end loop over predicted ATOF hits - } // end if - } // end loop over tracks + // RealMatrix measurementNoise = new Array2DRowRealMatrix( + // new double[][]{ + // {dr2, 0.0000, 0.0000}, + // {0.00, dphi2, 0.0000}, + // {0.00, 0.0000, dz2} + // });//3x3; + // hit.setMeasurementNoise(measurementNoise); + // ATOF_hits.put(trackid, hit); + + // IsHitSelected = true; + // break; + // } + // } // end loop over ATOF bank rows + // } //end loop over predicted ATOF hits + // } // end if + // } // end loop over tracks /// Clean AHDC bad hits double sigma = 0.5; // mm From 658e9fe27f16e53154da4ec0a130cbed21dcef3b Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 15:45:05 +0100 Subject: [PATCH 23/35] fill kftrack before propagation to ATOF --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 56 +++++++++---------- 1 file changed, 28 insertions(+), 28 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 60b0af550d..884a88fd32 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 @@ -167,6 +167,27 @@ public void propagation(ArrayList tracks, DataEvent event, final double m hit.setResidual(PostFitPropagator.residual(hit)); } } + + // Fill track and hit bank + // TO DO : s and p_drift have to be checked to be sure they represent what we want + double s = PostFitPropagator.stepper.sTot; + double p_drift = PostFitPropagator.stepper.p(); + int sum_adc = 0; + double sum_residuals = 0; + double chi2 = 0; + for (Hit hit : AHDC_hits) { + sum_adc += hit.getADC(); + sum_residuals += hit.getResidual(); + chi2 += Math.pow(hit.getResidual(),2)/hit.MeasurementNoiseMatrix().getEntry(0,0); + } + track.set_sum_adc(sum_adc); + track.set_sum_residuals(sum_residuals); + track.set_chi2(chi2/(AHDC_hits.size()-3)); + track.set_p_drift(p_drift); + track.set_dEdx(sum_adc/s); + track.set_path(s); + track.set_n_hits(AHDC_hits.size()); + // Projection towards the ATOF surfaces // R1 : lower surface of an ATOF bar // R2 : upper surface of an ATOF bar = lower surface of an ATOF wedge @@ -177,9 +198,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m double R1 = Math.hypot(pR1.x(), pR1.y()); double R2 = Math.hypot(pR2.x(), pR2.y()); double R3 = Math.hypot(pR3.x(), pR3.y()); - - System.out.println("ATOF surface R2 : " + R2); - System.out.println("ATOF surface R3 : " + R3); { // From last AHDC hit to surface R1 RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); @@ -189,8 +207,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m posR1.translateXYZ(0, 0, atof_alignement); int[] idR1 = predict_bar(ATOFdet, posR1); System.out.println("ATOF surface R1 : " + R1); - System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR1.x(), posR1.y(), posR1.z()); - System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR1[0], idR1[1], idR1[2]); + System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f) --> R (%.2f)\n", posR1.x(), posR1.y(), posR1.z(), Math.hypot(posR1.x(), posR1.y())); + System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR1[0], idR1[1], idR1[2]); // From surface R1 to surface R2 RadialSurfaceKFHit hitR2 = new RadialSurfaceKFHit(R2); PostFitPropagator.predict(hitR2, true); @@ -199,8 +217,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m posR2.translateXYZ(0, 0, atof_alignement); int[] idR2 = predict_wedge(ATOFdet, posR2); System.out.println("ATOF surface R2 : " + R2); - System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR2.x(), posR2.y(), posR2.z()); - System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR2[0], idR2[1], idR2[2]); + System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f)\n", posR2.x(), posR2.y(), posR2.z()); + System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR2[0], idR2[1], idR2[2]); // From surface R2 to surface R3 RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); PostFitPropagator.predict(hitR3, true); @@ -209,29 +227,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m posR3.translateXYZ(0, 0, atof_alignement); int[] idR3 = predict_wedge(ATOFdet, posR3); System.out.println("ATOF surface R3 : " + R3); - System.out.printf (" final position : x (.2f) y (.2f) z (.2f)\n", posR3.x(), posR3.y(), posR3.z()); - System.out.printf (" ---> sector (2d) layer (2d) component (2d)\n", idR3[0], idR3[1], idR3[2]); + System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f)\n", posR3.x(), posR3.y(), posR3.z()); + System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR3[0], idR3[1], idR3[2]); } - // Fill track and hit bank - // TO DO : s and p_drift have to be checked to be sure they represent what we want - double s = PostFitPropagator.stepper.sTot; - double p_drift = PostFitPropagator.stepper.p(); - int sum_adc = 0; - double sum_residuals = 0; - double chi2 = 0; - for (Hit hit : AHDC_hits) { - sum_adc += hit.getADC(); - sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2)/hit.MeasurementNoiseMatrix().getEntry(0,0); - } - track.set_sum_adc(sum_adc); - track.set_sum_residuals(sum_residuals); - track.set_chi2(chi2/(AHDC_hits.size()-3)); - track.set_p_drift(p_drift); - track.set_dEdx(sum_adc/s); - track.set_path(s); - track.set_n_hits(AHDC_hits.size()); + From 9f92d8f4da90215d39a3ef5bd0b32d8986a5c3cb Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 16:03:31 +0100 Subject: [PATCH 24/35] remove redundancy of Stepper in KFitter --- .../org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 9 ++++++--- .../jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 14 +++++++------- 2 files changed, 13 insertions(+), 10 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 03bf2894c4..84968345fc 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 @@ -21,7 +21,7 @@ public class KFitter { private RealVector stateEstimation; private RealMatrix errorCovariance; - public final Stepper stepper; + private Stepper stepper; private final Propagator propagator; private final HashMap materialHashMap; public double chi2 = 0; @@ -30,10 +30,10 @@ public class KFitter { private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; private double[] vertex_resolutions = {0.09,1e10}; // default values // dr^2 and dz^2 in mm^2 - public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator, final HashMap materialHashMap) { + public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Propagator propagator, final HashMap materialHashMap) { this.stateEstimation = initialStateEstimate; this.errorCovariance = initialErrorCovariance; - this.stepper = stepper; + this.stepper = new Stepper(initialStateEstimate.toArray()); this.propagator = propagator; this.materialHashMap = materialHashMap; } @@ -178,4 +178,7 @@ public void setVertexResolution(double[] res) { vertex_resolutions[1] = res[1]; } + /** Return a copy of the stepper. It is like a snapshot of the propagation. */ + public Stepper getStepper() {return new Stepper(stepper);} + } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index 884a88fd32..cc13ccae54 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 @@ -100,16 +100,15 @@ public void propagation(ArrayList tracks, DataEvent event, final double m ArrayList AHDC_hits = track.getHits(); Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java - // Start propagation - Stepper stepper = new Stepper(y); + // Initialize propagator RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B); Propagator propagator = new Propagator(RK4); // Initialization of the Kalman Fitter // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 - RealVector initialStateEstimate = new ArrayRealVector(stepper.y); + RealVector initialStateEstimate = new ArrayRealVector(y); RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); - KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); + KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, propagator, materialHashMap); if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); // Loop over number of iterations @@ -159,7 +158,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); // Post fit propagation (no correction) to set the residuals - KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); + KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Propagator(RK4), materialHashMap); // Projection towards AHDC hits for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); @@ -170,8 +169,9 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // Fill track and hit bank // TO DO : s and p_drift have to be checked to be sure they represent what we want - double s = PostFitPropagator.stepper.sTot; - double p_drift = PostFitPropagator.stepper.p(); + Stepper current_stepper = PostFitPropagator.getStepper(); + double s = current_stepper.sTot; + double p_drift = current_stepper.p(); int sum_adc = 0; double sum_residuals = 0; double chi2 = 0; From f855a2d1549193f1389f10021e79eb19a8b45371 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 16:23:20 +0100 Subject: [PATCH 25/35] define the hit of the beam once --- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 5 ---- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 25 ++++++++++--------- 2 files changed, 13 insertions(+), 17 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 84968345fc..e4ae341a21 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 @@ -172,11 +172,6 @@ public RealVector getStateEstimationVector() { public RealMatrix getErrorCovarianceMatrix() { return errorCovariance.copy(); } - - public void setVertexResolution(double[] res) { - vertex_resolutions[0] = res[0]; - vertex_resolutions[1] = res[1]; - } /** Return a copy of the stepper. It is like a snapshot of the propagation. */ public Stepper getStepper() {return new Stepper(stepper);} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index cc13ccae54..8f8e8d8df6 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 @@ -84,7 +84,17 @@ public void propagation(ArrayList tracks, DataEvent event, final double m row++; } } - + + // Define the beamline as an axtra hit for the track fitting + RadialKFHit hit_beam = new RadialKFHit(0, 0, vz_constraint); + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {vertex_resolutions[0], 0.0000, 0.0000}, + {0.00, 1e10, 0.0000}, + {0.00, 0.0000, vertex_resolutions[1]} + });//3x3; + hit_beam.setMeasurementNoise(measurementNoise); + // Loop over tracks for (Track track : tracks) { // Initialize state vector @@ -109,7 +119,6 @@ public void propagation(ArrayList tracks, DataEvent event, final double m RealVector initialStateEstimate = new ArrayRealVector(y); RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, propagator, materialHashMap); - if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); // Loop over number of iterations for (int k = 0; k < Niter; k++) { @@ -139,16 +148,8 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } // Backward propagation (first layer to beamline) { - RadialKFHit hit = new RadialKFHit(0, 0, vz_constraint); - RealMatrix measurementNoise = new Array2DRowRealMatrix( - new double[][]{ - {vertex_resolutions[0], 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, vertex_resolutions[1]} - });//3x3; - hit.setMeasurementNoise(measurementNoise); - TrackFitter.predict(hit, false); - TrackFitter.correct(hit); + TrackFitter.predict(hit_beam, false); + TrackFitter.correct(hit_beam); } } From b033b779faa510f2a4095018015cd7f311812cf0 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 16:50:32 +0100 Subject: [PATCH 26/35] clean --- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 1 - .../main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 1 - 2 files changed, 2 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 e4ae341a21..76b7f19f7d 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 @@ -28,7 +28,6 @@ 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 double[] vertex_resolutions = {0.09,1e10}; // default values // dr^2 and dz^2 in mm^2 public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Propagator propagator, final HashMap materialHashMap) { this.stateEstimation = initialStateEstimate; 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 8f8e8d8df6..1e1b499f81 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 @@ -21,7 +21,6 @@ import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; import org.jlab.geom.prim.Point3D; -import java.util.Random; /** From 88938f00c82dc268cc8dd9fdfbc8be52edfdbfa3 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 25 Feb 2026 19:02:27 +0100 Subject: [PATCH 27/35] keep the precision of the stepper at the end of the propagation --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java index b2b3818cb9..6f4f58a835 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java @@ -92,9 +92,8 @@ void propagate(Stepper stepper, KFHit hit, HashMap materialHas } else { // the distance between the step is not so big but the distance with respect to the hit starts to increase - // so, go back to the previous step (the best we have) and stop the propagation, set the default step size + // so, go back to the previous step (the best we have) and stop the propagation stepper.copyContent(prev_stepper); - stepper.h = 0.5; break; } } From b6bd41808abf9872e737da415db7a9f929699602 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 26 Feb 2026 14:12:19 +0100 Subject: [PATCH 28/35] add position, momentum, path of the track when crossing atof material --- etc/bankdefs/hipo4/alert.json | 68 +++++++++++++++++++ .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 26 +++++++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 39 ++++++++--- .../java/org/jlab/rec/ahdc/Track/Track.java | 26 +++++++ 4 files changed, 148 insertions(+), 11 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 64d4faeec1..afeaf6415a 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -408,6 +408,74 @@ "name": "sum_residuals", "type": "F", "info": "Sum of residuals (mm)" + }, { + "name": "atof_region", + "type": "I", + "info": "a flag to know if a track reach S1, S2, or S3; is n if Sn is reached; is 0 otherwise" + }, { + "name": "atof_s1_x", + "type": "F", + "info": "x position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_y", + "type": "F", + "info": "y position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_z", + "type": "F", + "info": "z position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_path", + "type": "F", + "info": "total path from the beamline to the surface s1" + }, { + "name": "atof_s1_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s1" + }, { + "name": "atof_s2_x", + "type": "F", + "info": "x position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_y", + "type": "F", + "info": "y position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_z", + "type": "F", + "info": "z position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_path", + "type": "F", + "info": "total path from the beamline to the surface s2" + }, { + "name": "atof_s2_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s2" + }, { + "name": "atof_s3_x", + "type": "F", + "info": "x position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_y", + "type": "F", + "info": "y position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_z", + "type": "F", + "info": "z position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_path", + "type": "F", + "info": "total path from the beamline to the surface s3" + }, { + "name": "atof_s3_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s3" + }, { + "name": "atof_component", + "type": "I", + "info": "sector*10000 + layer*100 + component : of the lower surface of a bar or the lower surface of a wedge depending or not the track reached these surfaces, priority is given to the wedges" } ] }, { 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 0884ab547f..c069a96ce0 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java @@ -7,6 +7,7 @@ import org.jlab.rec.ahdc.Cluster.Cluster; import org.jlab.rec.ahdc.DocaCluster.DocaCluster; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.Track.Track; @@ -153,6 +154,31 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { bank.setFloat("chi2", row, (float) track.get_chi2()); bank.setFloat("sum_residuals", row, (float) track.get_sum_residuals()); + // track projection on ATOF surface S1 + Stepper stepper_s1 = track.get_ATOF_S1_stepper(); + bank.setFloat("atof_s1_x", row, (float) stepper_s1.y[0]); + bank.setFloat("atof_s1_y", row, (float) stepper_s1.y[1]); + bank.setFloat("atof_s1_z", row, (float) stepper_s1.y[2]); + bank.setFloat("atof_s1_path", row, (float) stepper_s1.sTot); + bank.setFloat("atof_s1_p", row, (float) stepper_s1.p()); + // track projection on ATOF surface S2 + Stepper stepper_s2 = track.get_ATOF_S2_stepper(); + bank.setFloat("atof_s2_x", row, (float) stepper_s2.y[0]); + bank.setFloat("atof_s2_y", row, (float) stepper_s2.y[1]); + bank.setFloat("atof_s2_z", row, (float) stepper_s2.y[2]); + bank.setFloat("atof_s2_path", row, (float) stepper_s2.sTot); + bank.setFloat("atof_s2_p", row, (float) stepper_s2.p()); + // track projection on ATOF surface S3 + Stepper stepper_s3 = track.get_ATOF_S3_stepper(); + bank.setFloat("atof_s3_x", row, (float) stepper_s3.y[0]); + bank.setFloat("atof_s3_y", row, (float) stepper_s3.y[1]); + bank.setFloat("atof_s3_z", row, (float) stepper_s3.y[2]); + bank.setFloat("atof_s3_path", row, (float) stepper_s3.sTot); + bank.setFloat("atof_s3_p", row, (float) stepper_s3.p()); + + bank.setInt("atof_region", row, track.get_ATOF_region()); + bank.setInt("atof_component", row, track.get_ATOF_component()); + row++; } 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 1e1b499f81..e243610c74 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 @@ -47,7 +47,6 @@ public class KalmanFilter { private int counter = 0; // number of utilisation of the Kalman Filter HashMap ATOF_hits = null; HashMap> ATOF_hits_predicted = new HashMap<>(); - //ArrayList ATOF_hits_predicted = new ArrayList<>(); // list of sector, layer, component (for now : only wedges) AlertTOFDetector ATOFdet = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -202,34 +201,52 @@ public void propagation(ArrayList tracks, DataEvent event, final double m // From last AHDC hit to surface R1 RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); PostFitPropagator.predict(hitR1, true); + Stepper stepperR1 = PostFitPropagator.getStepper(); double[] vecR1 = PostFitPropagator.getStateEstimationVector().toArray(); Point3D posR1 = new Point3D(vecR1[0], vecR1[1], vecR1[2]); posR1.translateXYZ(0, 0, atof_alignement); int[] idR1 = predict_bar(ATOFdet, posR1); - System.out.println("ATOF surface R1 : " + R1); - System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f) --> R (%.2f)\n", posR1.x(), posR1.y(), posR1.z(), Math.hypot(posR1.x(), posR1.y())); - System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR1[0], idR1[1], idR1[2]); + track.set_ATOF_S1_stepper(stepperR1); + if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { + track.set_ATOF_region(1); + if (counter < 2) track.set_ATOF_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); + } + + + // System.out.println("###########################################################"); + // System.out.println("ATOF surface R1 : " + R1); + // System.out.printf (" final position : x (%f) y (%f) z (%f) --> R (%f)\n", posR1.x(), posR1.y(), posR1.z(), Math.hypot(posR1.x(), posR1.y())); + // System.out.println(" stepper.h : " + PostFitPropagator.getStepper().h); + // System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR1[0], idR1[1], idR1[2]); + // From surface R1 to surface R2 RadialSurfaceKFHit hitR2 = new RadialSurfaceKFHit(R2); PostFitPropagator.predict(hitR2, true); + Stepper stepperR2 = PostFitPropagator.getStepper(); double[] vecR2 = PostFitPropagator.getStateEstimationVector().toArray(); Point3D posR2 = new Point3D(vecR2[0], vecR2[1], vecR2[2]); posR2.translateXYZ(0, 0, atof_alignement); int[] idR2 = predict_wedge(ATOFdet, posR2); - System.out.println("ATOF surface R2 : " + R2); - System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f)\n", posR2.x(), posR2.y(), posR2.z()); - System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR2[0], idR2[1], idR2[2]); + track.set_ATOF_S2_stepper(stepperR2); + if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { + track.set_ATOF_region(2); + if (counter < 2) track.set_ATOF_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); + } // From surface R2 to surface R3 RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); PostFitPropagator.predict(hitR3, true); + Stepper stepperR3 = PostFitPropagator.getStepper(); double[] vecR3 = PostFitPropagator.getStateEstimationVector().toArray(); Point3D posR3 = new Point3D(vecR3[0], vecR3[1], vecR3[2]); posR3.translateXYZ(0, 0, atof_alignement); int[] idR3 = predict_wedge(ATOFdet, posR3); - System.out.println("ATOF surface R3 : " + R3); - System.out.printf (" final position : x (%.2f) y (%.2f) z (%.2f)\n", posR3.x(), posR3.y(), posR3.z()); - System.out.printf (" ---> sector (%2d) layer (%2d) component (%2d)\n", idR3[0], idR3[1], idR3[2]); - } + track.set_ATOF_S3_stepper(stepperR3); + if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { + track.set_ATOF_region(3); + //if (counter < 2) track.set_ATOF_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); + } + + } // end propagation towards ATOF surface 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 75a33d7397..5324ce8b37 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 @@ -7,6 +7,7 @@ import org.jlab.rec.ahdc.Cluster.Cluster; import org.jlab.rec.ahdc.HelixFit.HelixFitObject; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.PreCluster.PreClusterFinder; import org.jlab.rec.ahdc.AI.PreClustering; @@ -45,6 +46,19 @@ public class Track { {0.0 , 0.0 , 0.0 , 100 , 0.0 , 0.0}, {0.0 , 0.0 , 0.0 , 0.0 , 100 , 0.0}, {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 900}}); + // Position and momentum when the track crosses the ATOF surface + // S1 : lower surface of an ATOF bar + // S2 : upper surface of an ATOF bar = lower surface of an ATOF wedge + // S3 : upper surface of an ATOF wedge + Stepper ATOF_S1_stepper; + Stepper ATOF_S2_stepper; + Stepper ATOF_S3_stepper; + double ATOF_S1_radius; + double ATOF_S2_radius; + double ATOF_S3_radius; + int ATOF_region = 0; // is n if the trach reaches Sn, 0 otherwise (i.e does not reach S1) + int ATOF_component = 0; // predicted wedge (or bar if region is 1) just after the first fit (KF1) + // AHDC::aiprediction private int predicted_ATOF_sector = -1; @@ -210,4 +224,16 @@ public void set_trackId(int _trackId) { public int get_predicted_ATOF_layer() {return predicted_ATOF_layer;} public int get_predicted_ATOF_wedge() {return predicted_ATOF_wedge;} + // Projection of the Track on the ATOF surfaces + public void set_ATOF_S1_stepper(Stepper _stepper) {this.ATOF_S1_stepper = _stepper;} + public void set_ATOF_S2_stepper(Stepper _stepper) {this.ATOF_S2_stepper = _stepper;} + public void set_ATOF_S3_stepper(Stepper _stepper) {this.ATOF_S3_stepper = _stepper;} + public void set_ATOF_region(int _n) {this.ATOF_region = _n;} + public void set_ATOF_component(int _comp) {this.ATOF_component = _comp;} + public Stepper get_ATOF_S1_stepper() {return this.ATOF_S1_stepper;} + public Stepper get_ATOF_S2_stepper() {return this.ATOF_S2_stepper;} + public Stepper get_ATOF_S3_stepper() {return this.ATOF_S3_stepper;} + public int get_ATOF_region() {return this.ATOF_region;} + public int get_ATOF_component() {return this.ATOF_component;} + } From c3aaf8163e8b40579e3776d2b167e0ed5bd27e8a Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 27 Feb 2026 13:03:57 +0100 Subject: [PATCH 29/35] fix bugs --- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 38 +++++++++++-------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 19 +++++++--- .../org/jlab/service/alert/ALERTEngine.java | 30 ++++++++------- 3 files changed, 53 insertions(+), 34 deletions(-) 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 c069a96ce0..2557259f76 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 @@ -156,25 +156,33 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { // track projection on ATOF surface S1 Stepper stepper_s1 = track.get_ATOF_S1_stepper(); - bank.setFloat("atof_s1_x", row, (float) stepper_s1.y[0]); - bank.setFloat("atof_s1_y", row, (float) stepper_s1.y[1]); - bank.setFloat("atof_s1_z", row, (float) stepper_s1.y[2]); - bank.setFloat("atof_s1_path", row, (float) stepper_s1.sTot); - bank.setFloat("atof_s1_p", row, (float) stepper_s1.p()); + if (stepper_s1 != null) { + bank.setFloat("atof_s1_x", row, (float) stepper_s1.y[0]); + bank.setFloat("atof_s1_y", row, (float) stepper_s1.y[1]); + bank.setFloat("atof_s1_z", row, (float) stepper_s1.y[2]); + bank.setFloat("atof_s1_path", row, (float) stepper_s1.sTot); + bank.setFloat("atof_s1_p", row, (float) stepper_s1.p()); + } + // track projection on ATOF surface S2 Stepper stepper_s2 = track.get_ATOF_S2_stepper(); - bank.setFloat("atof_s2_x", row, (float) stepper_s2.y[0]); - bank.setFloat("atof_s2_y", row, (float) stepper_s2.y[1]); - bank.setFloat("atof_s2_z", row, (float) stepper_s2.y[2]); - bank.setFloat("atof_s2_path", row, (float) stepper_s2.sTot); - bank.setFloat("atof_s2_p", row, (float) stepper_s2.p()); + if (stepper_s2 != null) { + bank.setFloat("atof_s2_x", row, (float) stepper_s2.y[0]); + bank.setFloat("atof_s2_y", row, (float) stepper_s2.y[1]); + bank.setFloat("atof_s2_z", row, (float) stepper_s2.y[2]); + bank.setFloat("atof_s2_path", row, (float) stepper_s2.sTot); + bank.setFloat("atof_s2_p", row, (float) stepper_s2.p()); + } + // track projection on ATOF surface S3 Stepper stepper_s3 = track.get_ATOF_S3_stepper(); - bank.setFloat("atof_s3_x", row, (float) stepper_s3.y[0]); - bank.setFloat("atof_s3_y", row, (float) stepper_s3.y[1]); - bank.setFloat("atof_s3_z", row, (float) stepper_s3.y[2]); - bank.setFloat("atof_s3_path", row, (float) stepper_s3.sTot); - bank.setFloat("atof_s3_p", row, (float) stepper_s3.p()); + if (stepper_s3 != null) { + bank.setFloat("atof_s3_x", row, (float) stepper_s3.y[0]); + bank.setFloat("atof_s3_y", row, (float) stepper_s3.y[1]); + bank.setFloat("atof_s3_z", row, (float) stepper_s3.y[2]); + bank.setFloat("atof_s3_path", row, (float) stepper_s3.sTot); + bank.setFloat("atof_s3_p", row, (float) stepper_s3.p()); + } bank.setInt("atof_region", row, track.get_ATOF_region()); bank.setInt("atof_component", row, track.get_ATOF_component()); 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 e243610c74..df722d85fa 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 @@ -209,7 +209,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S1_stepper(stepperR1); if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { track.set_ATOF_region(1); - if (counter < 2) track.set_ATOF_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); + if (counter < 2 && idR1 != null) track.set_ATOF_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); } @@ -230,7 +230,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S2_stepper(stepperR2); if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { track.set_ATOF_region(2); - if (counter < 2) track.set_ATOF_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); + if (counter < 2 && idR2 != null) track.set_ATOF_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); } // From surface R2 to surface R3 RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); @@ -243,7 +243,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S3_stepper(stepperR3); if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { track.set_ATOF_region(3); - //if (counter < 2) track.set_ATOF_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); + //if (counter < 2 && idR3 != null) track.set_ATOF_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); } } // end propagation towards ATOF surface @@ -342,6 +342,7 @@ static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { wedge = c; } } + if (wedge == -1) return null; // find sector and layer int sector = -1; int layer = -1; @@ -356,7 +357,11 @@ static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { } } } - return new int[] {sector, layer, wedge}; + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, wedge}; + } } /** @@ -378,7 +383,11 @@ static public int[] predict_bar(AlertTOFDetector atof, Point3D pt) { } } } - return new int[] {sector, layer, 10}; + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, 10}; + } } public ArrayList get_adjacent_wedges(int[] identifiers) { diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index d7896d60a7..c703507085 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -353,20 +353,22 @@ public boolean processDataEvent(DataEvent event) { AHDC_hits.add(hit); } } - Track track = new Track(AHDC_hits); - // Initialise the position and the momentum using the information of the AHDC::track - // position : mm - // momentum : MeV - double x = trackBank.getFloat("x", row); - double y = trackBank.getFloat("y", row); - double z = trackBank.getFloat("z", row); - double px = trackBank.getFloat("px", row); - double py = trackBank.getFloat("py", row); - double pz = trackBank.getFloat("pz", row); - double[] vec = {x, y, z, px, py, pz}; - track.setPositionAndMomentumVec(vec); - track.set_trackId(trackid); - AHDC_tracks.add(track); + if (AHDC_hits.size() > 0) { + Track track = new Track(AHDC_hits); + // Initialise the position and the momentum using the information of the AHDC::track + // position : mm + // momentum : MeV + double x = trackBank.getFloat("x", row); + double y = trackBank.getFloat("y", row); + double z = trackBank.getFloat("z", row); + double px = trackBank.getFloat("px", row); + double py = trackBank.getFloat("py", row); + double pz = trackBank.getFloat("pz", row); + double[] vec = {x, y, z, px, py, pz}; + track.setPositionAndMomentumVec(vec); + track.set_trackId(trackid); + AHDC_tracks.add(track); + } } /// Intialise the Kalman Filter double magfieldfactor = runBank.getFloat("solenoid", 0); From bfd2471b4f61aef4a7f9411151c4664a5fdacc90 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 27 Feb 2026 14:01:56 +0100 Subject: [PATCH 30/35] store predicted atof components for each surfaces --- etc/bankdefs/hipo4/alert.json | 12 ++++++++++++ .../org/jlab/rec/ahdc/Banks/RecoBankWriter.java | 5 +++++ .../jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 4 +++- .../main/java/org/jlab/rec/ahdc/Track/Track.java | 16 +++++++++++++++- 4 files changed, 35 insertions(+), 2 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index afeaf6415a..d68f302790 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -476,6 +476,18 @@ "name": "atof_component", "type": "I", "info": "sector*10000 + layer*100 + component : of the lower surface of a bar or the lower surface of a wedge depending or not the track reached these surfaces, priority is given to the wedges" + }, { + "name": "atof_s1_comp", + "type": "I", + "info": "sector*10000 + layer*100 + component : atof component on s1" + }, { + "name": "atof_s2_comp", + "type": "I", + "info": "sector*10000 + layer*100 + component : atof component on s1" + }, { + "name": "atof_s3_comp", + "type": "I", + "info": "sector*10000 + layer*100 + component : atof component on s1" } ] }, { 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 2557259f76..145ed0a107 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 @@ -187,6 +187,11 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { bank.setInt("atof_region", row, track.get_ATOF_region()); bank.setInt("atof_component", row, track.get_ATOF_component()); + // tmp + bank.setInt("atof_s1_comp", row, track.get_ATOF_S1_component()); + bank.setInt("atof_s2_comp", row, track.get_ATOF_S2_component()); + bank.setInt("atof_s3_comp", row, track.get_ATOF_S3_component()); + row++; } 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 df722d85fa..63c0cbde3e 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 @@ -210,6 +210,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { track.set_ATOF_region(1); if (counter < 2 && idR1 != null) track.set_ATOF_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); + if (counter < 2 && idR1 != null) track.set_ATOF_S1_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); } @@ -231,6 +232,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { track.set_ATOF_region(2); if (counter < 2 && idR2 != null) track.set_ATOF_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); + if (counter < 2 && idR2 != null) track.set_ATOF_S2_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); } // From surface R2 to surface R3 RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); @@ -243,7 +245,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S3_stepper(stepperR3); if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { track.set_ATOF_region(3); - //if (counter < 2 && idR3 != null) track.set_ATOF_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); + if (counter < 2 && idR3 != null) track.set_ATOF_S3_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); } } // end propagation towards ATOF surface 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 5324ce8b37..566444b813 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 @@ -57,7 +57,12 @@ public class Track { double ATOF_S2_radius; double ATOF_S3_radius; int ATOF_region = 0; // is n if the trach reaches Sn, 0 otherwise (i.e does not reach S1) - int ATOF_component = 0; // predicted wedge (or bar if region is 1) just after the first fit (KF1) + int ATOF_component = -1; // predicted wedge (or bar if region is 1) just after the first fit (KF1) + + // tmp + int ATOF_S1_component = -1; + int ATOF_S2_component = -1; + int ATOF_S3_component = -1; // AHDC::aiprediction @@ -236,4 +241,13 @@ public void set_trackId(int _trackId) { public int get_ATOF_region() {return this.ATOF_region;} public int get_ATOF_component() {return this.ATOF_component;} + // tmp + public int get_ATOF_S1_component() {return this.ATOF_S1_component;} + public int get_ATOF_S2_component() {return this.ATOF_S2_component;} + public int get_ATOF_S3_component() {return this.ATOF_S3_component;} + public void set_ATOF_S1_component(int _comp) {this.ATOF_S1_component = _comp;} + public void set_ATOF_S2_component(int _comp) {this.ATOF_S2_component = _comp;} + public void set_ATOF_S3_component(int _comp) {this.ATOF_S3_component = _comp;} + + } From 34dcc689e2ea7dfb4f1d711b01de641739ccce77 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 27 Feb 2026 15:03:38 +0100 Subject: [PATCH 31/35] fix array access bug --- .../java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 8 +++++--- 1 file changed, 5 insertions(+), 3 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 63c0cbde3e..91f979eea4 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 @@ -133,9 +133,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m TrackFitter.predict(hit, true); TrackFitter.correct(hit); // Backward propagation to the last ahdc layer - Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); - TrackFitter.predict(hhit, false); - TrackFitter.correct(hhit); + if (AHDC_hits.size() > 0) { + Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); + TrackFitter.predict(hhit, false); + TrackFitter.correct(hhit); + } } } // Backward propagation (last layer to first layer) From eaa8f0adde26551843c4261200728bd2e857068e Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 2 Mar 2026 13:34:42 +0100 Subject: [PATCH 32/35] do not modify initial state and error matrix --- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 4 ++-- 1 file changed, 2 insertions(+), 2 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 76b7f19f7d..8fb0ecaafd 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 @@ -30,8 +30,8 @@ public class KFitter { private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Propagator propagator, final HashMap materialHashMap) { - this.stateEstimation = initialStateEstimate; - this.errorCovariance = initialErrorCovariance; + this.stateEstimation = initialStateEstimate.copy(); + this.errorCovariance = initialErrorCovariance.copy(); this.stepper = new Stepper(initialStateEstimate.toArray()); this.propagator = propagator; this.materialHashMap = materialHashMap; From 75d1b46f2da07209a5a8bc73feeba6b20d300779 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 2 Mar 2026 14:00:15 +0100 Subject: [PATCH 33/35] save error in x,y,z for the atof projections --- etc/bankdefs/hipo4/alert.json | 40 +++++++++++++++++-- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 22 +++++++++- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 21 ++++++---- .../java/org/jlab/rec/ahdc/Track/Track.java | 17 ++++++-- 4 files changed, 84 insertions(+), 16 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index d68f302790..99032b8c64 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -472,10 +472,6 @@ "name": "atof_s3_p", "type": "F", "info": "momentum whith which the track reaches the surface s3" - }, { - "name": "atof_component", - "type": "I", - "info": "sector*10000 + layer*100 + component : of the lower surface of a bar or the lower surface of a wedge depending or not the track reached these surfaces, priority is given to the wedges" }, { "name": "atof_s1_comp", "type": "I", @@ -488,6 +484,42 @@ "name": "atof_s3_comp", "type": "I", "info": "sector*10000 + layer*100 + component : atof component on s1" + }, { + "name": "atof_s1_sigma_x", + "type": "F", + "info": "error in atof_s1_x" + }, { + "name": "atof_s1_sigma_y", + "type": "F", + "info": "error in atof_s1_y" + }, { + "name": "atof_s1_sigma_z", + "type": "F", + "info": "error in atof_s1_z" + }, { + "name": "atof_s2_sigma_x", + "type": "F", + "info": "error in atof_s2_x" + }, { + "name": "atof_s2_sigma_y", + "type": "F", + "info": "error in atof_s2_y" + }, { + "name": "atof_s2_sigma_z", + "type": "F", + "info": "error in atof_s2_z" + }, { + "name": "atof_s3_sigma_x", + "type": "F", + "info": "error in atof_s3_x" + }, { + "name": "atof_s3_sigma_y", + "type": "F", + "info": "error in atof_s3_y" + }, { + "name": "atof_s3_sigma_z", + "type": "F", + "info": "error in atof_s3_z" } ] }, { 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 145ed0a107..675ecdb21e 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 @@ -1,5 +1,6 @@ package org.jlab.rec.ahdc.Banks; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.AI.InterCluster; @@ -185,13 +186,32 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { } bank.setInt("atof_region", row, track.get_ATOF_region()); - bank.setInt("atof_component", row, track.get_ATOF_component()); // tmp bank.setInt("atof_s1_comp", row, track.get_ATOF_S1_component()); bank.setInt("atof_s2_comp", row, track.get_ATOF_S2_component()); bank.setInt("atof_s3_comp", row, track.get_ATOF_S3_component()); + // error on atof_sn_(x,y,z) + RealMatrix atof_s1_errorMatrix = track.get_ATOF_S1_errorMatrix(); + RealMatrix atof_s2_errorMatrix = track.get_ATOF_S2_errorMatrix(); + RealMatrix atof_s3_errorMatrix = track.get_ATOF_S3_errorMatrix(); + if (atof_s1_errorMatrix != null) { + bank.setFloat("atof_s1_sigma_x", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(0, 0))); + bank.setFloat("atof_s1_sigma_y", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(1, 1))); + bank.setFloat("atof_s1_sigma_z", row, (float) Math.sqrt(atof_s1_errorMatrix.getEntry(2, 2))); + } + if (atof_s2_errorMatrix != null) { + bank.setFloat("atof_s2_sigma_x", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(0, 0))); + bank.setFloat("atof_s2_sigma_y", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(1, 1))); + bank.setFloat("atof_s2_sigma_z", row, (float) Math.sqrt(atof_s2_errorMatrix.getEntry(2, 2))); + } + if (atof_s3_errorMatrix != null) { + bank.setFloat("atof_s3_sigma_x", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(0, 0))); + bank.setFloat("atof_s3_sigma_y", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(1, 1))); + bank.setFloat("atof_s3_sigma_z", row, (float) Math.sqrt(atof_s3_errorMatrix.getEntry(2, 2))); + } + row++; } 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 91f979eea4..d06777748f 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 @@ -45,8 +45,8 @@ public class KalmanFilter { private double atof_alignement = -32.7; double vz_constraint = 0; private int counter = 0; // number of utilisation of the Kalman Filter - HashMap ATOF_hits = null; - HashMap> ATOF_hits_predicted = new HashMap<>(); + HashMap ATOF_hits = null; // trackid vs KFHit + HashMap> ATOF_hits_predicted = new HashMap<>(); // trackid vs (sector, layer, wedge) AlertTOFDetector ATOFdet = null; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -211,8 +211,10 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S1_stepper(stepperR1); if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { track.set_ATOF_region(1); - if (counter < 2 && idR1 != null) track.set_ATOF_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); - if (counter < 2 && idR1 != null) track.set_ATOF_S1_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); + if (counter < 2 && idR1 != null) { + track.set_ATOF_S1_component(10000*idR1[0] + 100*idR1[1] + idR1[2]); + track.set_ATOF_S1_errorMatrix(PostFitPropagator.getErrorCovarianceMatrix()); + } } @@ -233,8 +235,10 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S2_stepper(stepperR2); if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { track.set_ATOF_region(2); - if (counter < 2 && idR2 != null) track.set_ATOF_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); - if (counter < 2 && idR2 != null) track.set_ATOF_S2_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); + if (counter < 2 && idR2 != null) { + track.set_ATOF_S2_component(10000*idR2[0] + 100*idR2[1] + idR2[2]); + track.set_ATOF_S2_errorMatrix(PostFitPropagator.getErrorCovarianceMatrix()); + } } // From surface R2 to surface R3 RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); @@ -247,7 +251,10 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_ATOF_S3_stepper(stepperR3); if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { track.set_ATOF_region(3); - if (counter < 2 && idR3 != null) track.set_ATOF_S3_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); + if (counter < 2 && idR3 != null) { + track.set_ATOF_S3_component(10000*idR3[0] + 100*idR3[1] + idR3[2]); + track.set_ATOF_S3_errorMatrix(PostFitPropagator.getErrorCovarianceMatrix()); + } } } // end propagation towards ATOF surface 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 566444b813..3ec459ca4e 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 @@ -38,7 +38,7 @@ public class Track { private double dEdx = 0; ///< deposited energy per path length (adc/mm) private double p_drift = 0; ///< momentum in the drift region (MeV) private double path = 0; ///< length of the track (mm) - // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 + // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 (in the beamline) RealMatrix errorCovarianceMatrix = MatrixUtils.createRealMatrix(new double[][]{ {50 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0}, {0.0 , 50 , 0.0 , 0.0 , 0.0 , 0.0}, @@ -57,12 +57,16 @@ public class Track { double ATOF_S2_radius; double ATOF_S3_radius; int ATOF_region = 0; // is n if the trach reaches Sn, 0 otherwise (i.e does not reach S1) - int ATOF_component = -1; // predicted wedge (or bar if region is 1) just after the first fit (KF1) // tmp int ATOF_S1_component = -1; int ATOF_S2_component = -1; int ATOF_S3_component = -1; + + // tmp : compatible with stepper: need to access error in x,y,z + RealMatrix ATOF_S1_errorMatrix; + RealMatrix ATOF_S2_errorMatrix; + RealMatrix ATOF_S3_errorMatrix; // AHDC::aiprediction @@ -234,12 +238,10 @@ public void set_trackId(int _trackId) { public void set_ATOF_S2_stepper(Stepper _stepper) {this.ATOF_S2_stepper = _stepper;} public void set_ATOF_S3_stepper(Stepper _stepper) {this.ATOF_S3_stepper = _stepper;} public void set_ATOF_region(int _n) {this.ATOF_region = _n;} - public void set_ATOF_component(int _comp) {this.ATOF_component = _comp;} public Stepper get_ATOF_S1_stepper() {return this.ATOF_S1_stepper;} public Stepper get_ATOF_S2_stepper() {return this.ATOF_S2_stepper;} public Stepper get_ATOF_S3_stepper() {return this.ATOF_S3_stepper;} public int get_ATOF_region() {return this.ATOF_region;} - public int get_ATOF_component() {return this.ATOF_component;} // tmp public int get_ATOF_S1_component() {return this.ATOF_S1_component;} @@ -249,5 +251,12 @@ public void set_trackId(int _trackId) { public void set_ATOF_S2_component(int _comp) {this.ATOF_S2_component = _comp;} public void set_ATOF_S3_component(int _comp) {this.ATOF_S3_component = _comp;} + public RealMatrix get_ATOF_S1_errorMatrix() {return this.ATOF_S1_errorMatrix;} + public RealMatrix get_ATOF_S2_errorMatrix() {return this.ATOF_S2_errorMatrix;} + public RealMatrix get_ATOF_S3_errorMatrix() {return this.ATOF_S3_errorMatrix;} + public void set_ATOF_S1_errorMatrix(RealMatrix _matrix) { this.ATOF_S1_errorMatrix = _matrix;} + public void set_ATOF_S2_errorMatrix(RealMatrix _matrix) { this.ATOF_S2_errorMatrix = _matrix;} + public void set_ATOF_S3_errorMatrix(RealMatrix _matrix) { this.ATOF_S3_errorMatrix = _matrix;} + } From 2239f75b3336f343211f2ea5bdf0c143ecad11a1 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 4 Mar 2026 11:54:11 +0100 Subject: [PATCH 34/35] remove getLine from the methods to be implemented --- .../alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java | 1 - .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 1 - .../java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java | 5 ----- .../org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java | 5 ----- 4 files changed, 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 d9f222aff4..5defe95d2f 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 @@ -108,7 +108,6 @@ public double getDoca() { return doca; } - @Override public Line3D getLine() { return wireLine; } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java index 6884e6f011..525273757b 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -10,7 +10,6 @@ * @author Felix Touchte Codjo */ public interface KFHit { - public Line3D getLine(); public double distance(Point3D point3D); public double getRadius(); /** Return the measurement encoded in this KFHit */ diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java index 7931cc51bd..220118255b 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -41,11 +41,6 @@ public RealMatrix MeasurementNoiseMatrix() { public void setMeasurementNoise(RealMatrix measurementNoise) { this.measurementNoise= measurementNoise; } - - @Override - public Line3D getLine() { - return line; - } @Override public double distance(Point3D point3D) { diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java index c3e026cec5..6e633877bb 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java @@ -33,11 +33,6 @@ public RealVector MeasurementVector() { public RealMatrix MeasurementNoiseMatrix() { return new Array2DRowRealMatrix(new double[][]{{1e-8}}); } - - @Override - public Line3D getLine() { - return null; - } @Override public double distance(Point3D point3D) { From bbd9512672d44a625510d6c050b972e36615025d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 4 Mar 2026 19:13:56 +0100 Subject: [PATCH 35/35] start revision of the kalman filter --- .../org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 1 - .../rec/ahdc/KalmanFilter/KalmanFilter.java | 142 ++++++++---------- .../java/org/jlab/rec/ahdc/Track/Track.java | 16 +- .../org/jlab/service/alert/ALERTEngine.java | 96 +++++++++--- 4 files changed, 156 insertions(+), 99 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java index 525273757b..3a2e693192 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -2,7 +2,6 @@ import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; -import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; /** * An interface to unify the hits used the Kalman Filter (e.g AHDC hits, ATOF hits, beamline) 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 d06777748f..936a29bd12 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java @@ -1,18 +1,14 @@ package org.jlab.rec.ahdc.KalmanFilter; import java.util.ArrayList; -import java.util.Collections; import java.util.HashMap; import java.util.Random; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.io.base.DataBank; -import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.Track; @@ -36,20 +32,21 @@ public class KalmanFilter { public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;} + HashMap materialHashMap = MaterialMap.generateMaterials(); + private PDGParticle particle; private int Niter = 40; // number of iterations for the Kalman Filter private boolean IsVtxDefined = false; // implemented but not used yet - private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} - // mm, they the misalignement with respect to the AHDC - private double clas_alignement = +54; - private double atof_alignement = -32.7; - double vz_constraint = 0; + private double vz_constraint = 0; + // mm, they are the misalignement with respect to the AHDC: the are defined in ALERTEngine + private double clas_alignement = 0; + private double atof_alignement = 0; + private int counter = 0; // number of utilisation of the Kalman Filter - HashMap ATOF_hits = null; // trackid vs KFHit HashMap> ATOF_hits_predicted = new HashMap<>(); // trackid vs (sector, layer, wedge) - AlertTOFDetector ATOFdet = null; + AlertTOFDetector ATOFdet = null; // a reference to the ATOF geometry - public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { + public void propagation(ArrayList tracks, final double magfield, boolean IsMC) { try { counter++; @@ -58,44 +55,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; - HashMap materialHashMap = MaterialMap.generateMaterials(); - // Recover the vertex of the electron - if (event.hasBank("REC::Particle") && !IsVtxDefined) { // as we run the KF several times, !IsVtxDefined prevent to look for the vertex again - DataBank recBank = event.getBank("REC::Particle"); - int row = 0; - while ((!IsVtxDefined) && row < recBank.rows()) { - if (recBank.getInt("pid", row) == 11) { - IsVtxDefined = true; - vz_constraint = 10*recBank.getFloat("vz",row) + (IsMC ? 0 : clas_alignement); // mm - - // TO DO: compute electron resolution as function of p and theta - // the fine tuning will be done later - //double px = recBank.getFloat("px",row); - //double py = recBank.getFloat("py",row); - //double pz = recBank.getFloat("pz",row); - //double p = Math.sqrt(px*px+py*py+pz*pz); - //double theta = Math.acos(pz/p); - - vertex_resolutions[0] = 0.09; - vertex_resolutions[1] = 64; - } - row++; - } - } - - // Define the beamline as an axtra hit for the track fitting - RadialKFHit hit_beam = new RadialKFHit(0, 0, vz_constraint); - RealMatrix measurementNoise = new Array2DRowRealMatrix( - new double[][]{ - {vertex_resolutions[0], 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, vertex_resolutions[1]} - });//3x3; - hit_beam.setMeasurementNoise(measurementNoise); + // Loop over tracks for (Track track : tracks) { - // Initialize state vector + /// Initialize state vector double x0 = 0.0; double y0 = 0.0; double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); @@ -104,71 +68,82 @@ public void propagation(ArrayList tracks, DataEvent event, final double m double pz0 = track.get_pz(); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; - // Read list of hits + + /// Read list of hits + RadialKFHit beam_hit = track.getBeamlineHit(); ArrayList AHDC_hits = track.getHits(); - Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java + ArrayList ATOF_hits = track.getATOFHits(); - // Initialize propagator + /// Initialize propagator RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B); Propagator propagator = new Propagator(RK4); - // Initialization of the Kalman Fitter - // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 + /// Initialization of the Kalman Fitter RealVector initialStateEstimate = new ArrayRealVector(y); RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, propagator, materialHashMap); // Loop over number of iterations for (int k = 0; k < Niter; k++) { - // Forward propagation + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); } - // Forward propagation towards the ATOF bar - { - RadialKFHit hit = ATOF_hits.get(track.get_trackId()); - if (hit != null) { + + // Take into account ATOF hits + if (!AHDC_hits.isEmpty() && !ATOF_hits.isEmpty()) { + // Forward propagation in ATOF + for (KFHit hit : ATOF_hits) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); - // Backward propagation to the last ahdc layer - if (AHDC_hits.size() > 0) { - Hit hhit = AHDC_hits.get(AHDC_hits.size()-1); - TrackFitter.predict(hhit, false); - TrackFitter.correct(hhit); - } + } + // Backward propagation in ATOF + for (int i = ATOF_hits.size()-2; i >= 0; i--){ + KFHit hit = ATOF_hits.get(i); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); + } + // Backward propagation to the last AHDC layer + { + Hit hit = AHDC_hits.getLast(); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); } } - // Backward propagation (last layer to first layer) - for (int i = AHDC_hits.size() - 2; i >= 0; i--) { + + // Backward propagation (from AHDC last layer to first AHDC layer) + for (int i = AHDC_hits.size()-2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } - // Backward propagation (first layer to beamline) + + // Backward propagation (from first AHDC layer to beamline) { - TrackFitter.predict(hit_beam, false); - TrackFitter.correct(hit_beam); + TrackFitter.predict(beam_hit, false); + TrackFitter.correct(beam_hit); } - } + } // end loop over nb. iteration + /// Output: position and momentum RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); - // Post fit propagation (no correction) to set the residuals + /// Post fit propagation (no correction) KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Propagator(RK4), materialHashMap); - // Projection towards AHDC hits + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); - if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits - hit.setResidual(PostFitPropagator.residual(hit)); - } + hit.setResidual(PostFitPropagator.residual(hit)); // output: residual } - // Fill track and hit bank + // Fill values for AHDC::kftrack // TO DO : s and p_drift have to be checked to be sure they represent what we want Stepper current_stepper = PostFitPropagator.getStepper(); double s = current_stepper.sTot; @@ -189,6 +164,11 @@ public void propagation(ArrayList tracks, DataEvent event, final double m track.set_path(s); track.set_n_hits(AHDC_hits.size()); + + //////// + /// Revision of the code: I'm here + /// ///// + // Projection towards the ATOF surfaces // R1 : lower surface of an ATOF bar // R2 : upper surface of an ATOF bar = lower surface of an ATOF wedge @@ -199,7 +179,7 @@ public void propagation(ArrayList tracks, DataEvent event, final double m double R1 = Math.hypot(pR1.x(), pR1.y()); double R2 = Math.hypot(pR2.x(), pR2.y()); double R3 = Math.hypot(pR3.x(), pR3.y()); - { + { // From last AHDC hit to surface R1 RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); PostFitPropagator.predict(hitR1, true); @@ -270,11 +250,9 @@ public void propagation(ArrayList tracks, DataEvent event, final double m } } public void set_Niter(int Niter) {this.Niter = Niter;} - public int get_Niter() {return this.Niter;} - public void set_particle(PDGParticle particle) {this.particle = particle;} + public int get_Niter() {return this.Niter;} + public void set_particle(PDGParticle particle) {this.particle = particle;} public PDGParticle get_particle() {return this.particle;} - public void set_ATOF_hits(HashMap ATOF_hits){ this.ATOF_hits = ATOF_hits;}; - public HashMap get_ATOF_hits() { return this.ATOF_hits;} public HashMap> get_ATOF_hits_predicted() {return this.ATOF_hits_predicted;} // Test @@ -441,4 +419,8 @@ else if (layer == 3) { } public void set_ATOF_detector(AlertTOFDetector atof) { this.ATOFdet = atof;} + public void set_clas_alignement(double _shift) {this.clas_alignement = _shift;} + public void set_atof_alignement(double _shift) {this.atof_alignement = _shift;} + public void set_vz_constraint(double _vz) {this.vz_constraint = _vz;} + public void set_vertex_flag(boolean _flag) {this.IsVtxDefined = _flag;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java index 3ec459ca4e..969cb7eaab 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 @@ -7,6 +7,7 @@ import org.jlab.rec.ahdc.Cluster.Cluster; import org.jlab.rec.ahdc.HelixFit.HelixFitObject; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; import org.jlab.rec.ahdc.KalmanFilter.Stepper; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.PreCluster.PreClusterFinder; @@ -21,7 +22,9 @@ public class Track { private List _Clusters = new ArrayList<>(); private List _InterClusters = new ArrayList<>(); private boolean _Used = false; - private final ArrayList hits = new ArrayList<>(); + private final ArrayList hits = new ArrayList<>(); // AHDC hits + private ArrayList ATOF_hits = new ArrayList<>(); + private RadialKFHit beamline_hit = null; private int trackId = -1; ///< id of the track private int n_hits = 0; ///< number of hits @@ -142,6 +145,17 @@ public ArrayList getHits() { return hits; } + public ArrayList getATOFHits() { + return this.ATOF_hits; + } + + public RadialKFHit getBeamlineHit() { + return this.beamline_hit; + } + + public void setATOFHits(ArrayList _ATOF_hits) {this.ATOF_hits = _ATOF_hits;} + public void setBeamlineHit(RadialKFHit _beamline_hit) {this.beamline_hit = _beamline_hit;} + @Override public String toString() { return "Track{" + "_Clusters=" + _Clusters + '}'; diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index c703507085..4a560ae110 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -4,6 +4,7 @@ import ai.djl.translate.TranslateException; import java.io.File; import java.util.ArrayList; +import java.util.Collections; import java.util.Iterator; import java.util.concurrent.atomic.AtomicInteger; import java.util.HashMap; @@ -326,8 +327,42 @@ public boolean processDataEvent(DataEvent event) { /// Kalmam Filter /// /////////////////////////////////////// - /// Read the list of tracks/hits from the banks AHDC::track and AHDC::hits + /// Pre conditions if (!event.hasBank("AHDC::track")) {return false;} + if (!event.hasBank("AHDC::hits")) {return false;} + + /// tmp: misalignement with respect to the center of the AHDC (mm) + double clas_alignement = +54; + double atof_alignement = -32.7; + + /// Read the electron vertex + double vz_electron = 0; + double[] vz_error2 = {0.09, 1e10, 1e10}; // mm^2, radians^2, mm^2, error on r, phi, z + boolean IsVertexDefined = false; + if (event.hasBank("REC::Particle")) { + DataBank recBank = event.getBank("REC::Particle"); + for (int row = 0; row < recBank.rows(); row++) { + if (recBank.getInt("pid", row) == 11) { + vz_electron = 10*recBank.getFloat("vz",row); // conversion in mm + IsVertexDefined = true; + + + //double px = recBank.getFloat("px",row); + //double py = recBank.getFloat("py",row); + //double pz = recBank.getFloat("pz",row); + //double p = Math.sqrt(px*px+py*py+pz*pz); + //double theta = Math.acos(pz/p); + + // set the resolutions on r and z! to be done + vz_error2[0] = 0.09; // should depend on p and theta + vz_error2[2] = 64; // should depend on p and theta + + break; // only look at the first electron + } + } + } + + /// Read the list of tracks/hits from the banks AHDC::track and AHDC::hits DataBank trackBank = event.getBank("AHDC::track"); DataBank hitBank = event.getBank("AHDC::hits"); ArrayList AHDC_tracks = new ArrayList<>(); @@ -344,7 +379,7 @@ public boolean processDataEvent(DataEvent event) { double doca = hitBank.getDouble("doca", hit_row); double time = hitBank.getDouble("time", hit_row); double tot = hitBank.getDouble("timeOverThreshold", hit_row); - // warning : adc is the calibrated one, we need the adc for the Kalman filter + // warning : adc is the calibrated one, we need the adc for the Kalman filter ! Hit hit = new Hit(id, superlayer, layer, wire, doca, adc, time); hit.setWirePosition(AHDC); hit.setTrackId(trackid); @@ -354,6 +389,7 @@ public boolean processDataEvent(DataEvent event) { } } if (AHDC_hits.size() > 0) { + Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java Track track = new Track(AHDC_hits); // Initialise the position and the momentum using the information of the AHDC::track // position : mm @@ -370,16 +406,24 @@ public boolean processDataEvent(DataEvent event) { AHDC_tracks.add(track); } } - /// Intialise the Kalman Filter - double magfieldfactor = runBank.getFloat("solenoid", 0); - double magfield = 50*magfieldfactor; + + /// Associate the electron vertex (the beamline hit) to each track boolean IsMC = event.hasBank("MC::Particle"); - PDGParticle proton = PDGDatabase.getParticleById(2212); - int Niter = 40; - KalmanFilter KF = new KalmanFilter(proton, Niter); + double vz_constraint = vz_electron + (IsMC ? 0 : clas_alignement); // we don't have the misalignment in simulation + for (Track track : AHDC_tracks) { + RadialKFHit hit_beam = new RadialKFHit(0, 0, vz_constraint); + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {vz_error2[0], 0.0000 , 0.0000}, + {0.0000 , vz_error2[1], 0.0000}, + {0.0000 , 0.0000 , vz_error2[2]} + });//3x3; + hit_beam.setMeasurementNoise(measurementNoise); + track.setBeamlineHit(hit_beam); + } - /// Add ATOF wedge hits - HashMap ATOF_hits = new HashMap<>(); + /// Look for ATOF wedge hits predicted by the AI + HashMap map_ATOF_hits = new HashMap<>(); for (Pair pair : matched_ATOF_hit_id) { int trackid = pair.getKey(); int atofid = pair.getValue(); @@ -390,7 +434,7 @@ public boolean processDataEvent(DataEvent event) { double x = bank_ATOFHits.getFloat("x", row); double y = bank_ATOFHits.getFloat("y", row); double z = bank_ATOFHits.getFloat("z", row); - z -= 32.3; // there is a shift between AHDC and ATOF (still don't know why) ! + z += atof_alignement; // there is a shift between AHDC and ATOF (still don't know why) ! RadialKFHit hit = new RadialKFHit(x, y, z); // error on r double wedge_width = 20; //mm @@ -409,16 +453,34 @@ public boolean processDataEvent(DataEvent event) { {0.00, 0.0000, dz2} });//3x3; hit.setMeasurementNoise(measurementNoise); - ATOF_hits.put(trackid, hit); + map_ATOF_hits.put(trackid, hit); } } } } - KF.set_ATOF_hits(ATOF_hits); - KF.set_ATOF_detector(ATOF); - /// First propagation : each AHDC_tracks will be fitted - KF.propagation(AHDC_tracks, event, magfield, IsMC); + /// Associate the ATOF hits to each track + for (Track track : AHDC_tracks) { + RadialKFHit hit = map_ATOF_hits.get(track.get_trackId()); + ArrayList list = new ArrayList<>(); + if (hit != null) list.add(hit); // for now, we only consider one hit in the ATOF + track.setATOFHits(list); + } + + /// Intialise the Kalman Filter + double magfieldfactor = runBank.getFloat("solenoid", 0); + double magfield = 50*magfieldfactor; + PDGParticle proton = PDGDatabase.getParticleById(2212); + int Niter = 40; + KalmanFilter KF = new KalmanFilter(proton, Niter); + KF.set_ATOF_detector(ATOF); // Reference the ATOF geometry in the Kalman Filter + KF.set_clas_alignement(clas_alignement); + KF.set_atof_alignement(atof_alignement); + KF.set_vz_constraint(vz_constraint); + KF.set_vertex_flag(IsVertexDefined); + + /// Do a first propagation + KF.propagation(AHDC_tracks, magfield, IsMC); // // Get ATOF wedges predicted // HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); @@ -483,7 +545,7 @@ public boolean processDataEvent(DataEvent event) { /// Second propagation : each AHDC_tracks will be fitted KF.set_Niter(15); - KF.propagation(AHDC_tracks, event, magfield, IsMC); + KF.propagation(AHDC_tracks, magfield, IsMC); /// write the AHDC::kftrack bank in the event org.jlab.rec.ahdc.Banks.RecoBankWriter ahdc_writer = new org.jlab.rec.ahdc.Banks.RecoBankWriter();