From 1f21a44c8d0a13612f8302f793e6dae740b9791e Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 4 Nov 2025 11:01:21 -0500 Subject: [PATCH 01/34] use the same time2distance for simu and real data We now invert the ccdb time2distance in simulation. --- .../alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java | 4 ---- 1 file changed, 4 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java index 54a8c047a6..8bdddfeff7 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java @@ -80,10 +80,6 @@ public final void fetch_AHDCHits(DataEvent event, AlertDCDetector detector) { // Hit selection : wfType and additional cuts if (((wfType <= 2) && (adc >= adc_min) && (adc <= adc_max) && (time >= t_min) && (time <= t_max) && (timeOverThreshold >= tot_min) && (timeOverThreshold <= tot_max) && (adcOffset >= ped_min) && (adcOffset <= ped_max)) || sim) { double doca = p0 + p1*Math.pow(time,1.0) + p2*Math.pow(time,2.0) + p3*Math.pow(time,3.0) + p4*Math.pow(time,4.0) + p5*Math.pow(time, 5.0); - if (sim) { - time += 5; // correction from mctime - systematic error of the decoding - doca = -0.0497 - 0.00667*Math.pow(time,1.0) + 0.389*Math.pow(time,1.0/2) - 0.189*Math.pow(time,1.0/3); - } if (time < 0) doca = 0; Hit h = new Hit(id, superlayer, layer, wire, doca, adc, time); h.setWirePosition(detector); From 5818e7c78456be4523b84d34d6faab1e574cca82 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 6 Nov 2025 09:36:43 -0500 Subject: [PATCH 02/34] define KFMonitor --- .../org/jlab/rec/ahdc/Track/KFMonitor.java | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java new file mode 100644 index 0000000000..3cceafb856 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java @@ -0,0 +1,26 @@ +package org.jlab.rec.ahdc.Track; + +import org.apache.commons.math3.linear.RealVector; + +public class KFMonitor { + private int Niter; //< iteration number of the Kalman Filter algorithm + private int orientation; //< forward (0), backward (1), postfit (2) propagation + private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 + private int status; //< state just after: prediction (0) or correction (1) + private RealVector state; //< vector containing x, y, z, px, py, pz + + // constructor + public KFMonitor(int _Niter, int _orientation, int _indicator, int _status, RealVector _state) { + this.Niter = _Niter; + this.orientation = _orientation; + this.indicator = _indicator; + this.status = _status; + this.state = _state.copy(); + } + + public int get_Niter() {return Niter;} + public int get_orientation() {return orientation;} + public int get_indicator() {return indicator;} + public int get_status() {return status;} + public RealVector get_state() {return state.copy();} +} From d7ecc36030965d54be237dd46e2f96b57aaf8fe3 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 6 Nov 2025 10:24:41 -0500 Subject: [PATCH 03/34] store kalman filter steps in a list of KFMonitor --- .../jlab/rec/ahdc/KalmanFilter/Indicator.java | 24 +++++++++++++++++++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 6 +++++ .../org/jlab/rec/ahdc/Track/KFMonitor.java | 2 +- .../java/org/jlab/rec/ahdc/Track/Track.java | 7 ++++++ 4 files changed, 38 insertions(+), 1 deletion(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java index 58e23ec9ba..09f93054fe 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java @@ -23,4 +23,28 @@ public boolean haveAHit() { if (this.R == 0 && !direction) res = true; return res; } + + + // This method is not self consistent. We should not have hardcoded values + // It is just for testing purpose + public int getUniqueId() { + if (hit != null) { + return (hit.getSuperLayer()*10 + hit.getLayer())*100 + hit.getWire(); + } + else if (Math.abs(R) < 1e-9) { + // beamline + return 0; + } + else if (Math.abs(R - 3.0) < 1e-9) { + // inner face of the target straw + return 1; + } + else if (Math.abs(R - 3.060) < 1e-9) { + // outer face of the target straw + return 2; + } + else { + return -1; + } + } } 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 7c53a3397b..f186798e40 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 @@ -12,6 +12,7 @@ import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Track.Track; +import org.jlab.rec.ahdc.Track.KFMonitor; import java.util.ArrayList; import java.util.HashMap; @@ -168,6 +169,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double //TrackFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { TrackFitter.predict(indicator); + track.add_KFMonitor(new KFMonitor(k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { if( k==0 && indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ @@ -175,14 +177,17 @@ private void propagation(ArrayList tracks, DataEvent event, final double } } TrackFitter.correct(indicator); + track.add_KFMonitor(new KFMonitor(k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); } } //System.out.println("--------- BackWard propagation !! ---------"); for (Indicator indicator : backwardIndicators) { TrackFitter.predict(indicator); + track.add_KFMonitor(new KFMonitor(k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { TrackFitter.correct(indicator); + track.add_KFMonitor(new KFMonitor(k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); } } } @@ -195,6 +200,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); for (Indicator indicator : forwardIndicators) { PostFitPropagator.predict(indicator); + track.add_KFMonitor(new KFMonitor(Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { if( indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java index 3cceafb856..6223132bfd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java @@ -5,7 +5,7 @@ public class KFMonitor { private int Niter; //< iteration number of the Kalman Filter algorithm private int orientation; //< forward (0), backward (1), postfit (2) propagation - private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 + private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]] private int status; //< state just after: prediction (0) or correction (1) private RealVector state; //< vector containing x, y, z, px, py, pz 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 ea1548a083..279ccf0fd9 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 @@ -41,6 +41,9 @@ public class Track { private double dEdx_kf = 0; ///< deposited energy per path length (adc/mm) private double p_drift_kf = 0; ///< momentum in the drift region (MeV) private double path_kf = 0; ///< length of the track (mm) + // Monitoring of the Kalman Filter + private List ListOfKFMonitors = new ArrayList<>(); + public Track(List clusters) { this._Clusters = clusters; @@ -201,4 +204,8 @@ public double getPz0_kf() { public double get_p_drift_kf() {return p_drift_kf;} public double get_path_kf() {return path_kf;} + // KF monitoring + public void add_KFMonitor(KFMonitor _monitor) { ListOfKFMonitors.add(_monitor);} + public List get_ListOfKFMonitors() { return ListOfKFMonitors;} + } From 6b0350b7ba2a0c2f358eb2b69638285324a50e7d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 6 Nov 2025 10:39:29 -0500 Subject: [PATCH 04/34] add a new bank structure for AHDC KF track monitoring --- etc/bankdefs/hipo4/alert.json | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index dd57254c5c..e404dfc734 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -414,5 +414,23 @@ {"name": "y5", "type": "F", "info": "Y5 position of the 5th superprecluster (mm)"}, {"name": "pred", "type": "F", "info": "Prediction of the model: 0 mean bad track; 1 mean good track"} ] + }, + { + "name": "AHDC::kftrack:mon", + "group": 23000, + "item": 32, + "info": "Monitoring of the Kalman Filter a given kftrack", + "entries": [ + {"name": "Niter", "type": "S", "info": "iteration number of the Kalman Filter algorithm"}, + {"name": "orientation", "type": "S", "info": "forward (0), backward (1), postfit (2) propagation"}, + {"name": "indicator", "type": "S", "info": "wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]]"}, + {"name": "status", "type": "S", "info": "state just after: prediction (0) or correction (1)"}, + {"name": "x", "type": "F", "info": "x position of the state (mm)"}, + {"name": "y", "type": "F", "info": "y position of the state (mm)"}, + {"name": "z", "type": "F", "info": "z position of the state (mm)"}, + {"name": "px", "type": "F", "info": "px momentum of the state (MeV)"}, + {"name": "py", "type": "F", "info": "py momentum of the state (MeV)"}, + {"name": "pz", "type": "F", "info": "pz momentum of the state (MeV)"} + ] } ] From b8d421325b8118236e044c91feca73ccc0fc1808 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 6 Nov 2025 10:47:21 -0500 Subject: [PATCH 05/34] add trackid in KFMonitor --- etc/bankdefs/hipo4/alert.json | 1 + .../org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 10 +++++----- .../main/java/org/jlab/rec/ahdc/Track/KFMonitor.java | 5 ++++- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index e404dfc734..f7aa7f69b3 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -421,6 +421,7 @@ "item": 32, "info": "Monitoring of the Kalman Filter a given kftrack", "entries": [ + {"name": "trackid", "type": "S", "info": "trackid"}, {"name": "Niter", "type": "S", "info": "iteration number of the Kalman Filter algorithm"}, {"name": "orientation", "type": "S", "info": "forward (0), backward (1), postfit (2) propagation"}, {"name": "indicator", "type": "S", "info": "wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]]"}, 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 f186798e40..cc768f0112 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 @@ -169,7 +169,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double //TrackFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { TrackFitter.predict(indicator); - track.add_KFMonitor(new KFMonitor(k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { if( k==0 && indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ @@ -177,17 +177,17 @@ private void propagation(ArrayList tracks, DataEvent event, final double } } TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); } } //System.out.println("--------- BackWard propagation !! ---------"); for (Indicator indicator : backwardIndicators) { TrackFitter.predict(indicator); - track.add_KFMonitor(new KFMonitor(k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); } } } @@ -200,7 +200,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); for (Indicator indicator : forwardIndicators) { PostFitPropagator.predict(indicator); - track.add_KFMonitor(new KFMonitor(Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); if (indicator.haveAHit()) { if( indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java index 6223132bfd..190f9cbf1d 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java @@ -3,6 +3,7 @@ import org.apache.commons.math3.linear.RealVector; public class KFMonitor { + private int trackid; //< trackid private int Niter; //< iteration number of the Kalman Filter algorithm private int orientation; //< forward (0), backward (1), postfit (2) propagation private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]] @@ -10,7 +11,8 @@ public class KFMonitor { private RealVector state; //< vector containing x, y, z, px, py, pz // constructor - public KFMonitor(int _Niter, int _orientation, int _indicator, int _status, RealVector _state) { + public KFMonitor(int _trackid, int _Niter, int _orientation, int _indicator, int _status, RealVector _state) { + this.trackid = _trackid; this.Niter = _Niter; this.orientation = _orientation; this.indicator = _indicator; @@ -18,6 +20,7 @@ public KFMonitor(int _Niter, int _orientation, int _indicator, int _status, Real this.state = _state.copy(); } + public int get_trackid() {return trackid;} public int get_Niter() {return Niter;} public int get_orientation() {return orientation;} public int get_indicator() {return indicator;} From cc8aa0e08d4e9e1384e1e2624c47e436834c8c8a Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 6 Nov 2025 11:21:37 -0500 Subject: [PATCH 06/34] add AHDC::kftrack:mon in AHDCEngine --- etc/bankdefs/util/bankSplit.py | 2 +- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 37 +++++++++++++++++++ .../org/jlab/service/ahdc/AHDCEngine.java | 4 +- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/etc/bankdefs/util/bankSplit.py b/etc/bankdefs/util/bankSplit.py index 801902b9ca..4ae0f4dfa9 100755 --- a/etc/bankdefs/util/bankSplit.py +++ b/etc/bankdefs/util/bankSplit.py @@ -65,7 +65,7 @@ def create(dirname, banklist): raster = ["RASTER::position"] rich = ["RICH::tdc","RICH::Ring","RICH::Particle"] rtpc = ["RTPC::hits","RTPC::tracks","RTPC::KFtracks"] -alert = ["AHDC::track", "AHDC::mc", "AHDC::hits", "AHDC::preclusters", "AHDC::clusters", "AHDC::kftrack", "AHDC::ai:prediction"] +alert = ["AHDC::track", "AHDC::mc", "AHDC::hits", "AHDC::preclusters", "AHDC::clusters", "AHDC::kftrack", "AHDC::ai:prediction", "AHDC::kftrack:mon"] dets = band + raster + rich + rtpc + alert # additions for the calibration schema: 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 38a701cd3b..b6d9c37e6a 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,8 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.Track.Track; +import org.jlab.rec.ahdc.Track.KFMonitor; +import org.apache.commons.math3.linear.RealVector; import java.util.ArrayList; @@ -179,4 +181,39 @@ public DataBank fillAIPrediction(DataEvent event, ArrayList pre return bank; } + + public DataBank fillKFMonitorBank(DataEvent event, ArrayList tracks) { + + int nentries = 0; + for (Track track : tracks) { + nentries += track.get_ListOfKFMonitors().size(); + } + + DataBank bank = event.createBank("AHDC::kftrack:mon", nentries); + + int row = 0; + + for (Track track : tracks) { + if (track == null) continue; + for (KFMonitor monitor : track.get_ListOfKFMonitors()) { + RealVector state = monitor.get_state(); + bank.setFloat("x", row, (float) state.getEntry(0)); + bank.setFloat("y", row, (float) state.getEntry(1)); + bank.setFloat("z", row, (float) state.getEntry(2)); + bank.setFloat("px", row, (float) state.getEntry(3)); + bank.setFloat("py", row, (float) state.getEntry(4)); + bank.setFloat("pz", row, (float) state.getEntry(5)); + bank.setShort("trackid", row, (short) monitor.get_trackid()); + bank.setShort("Niter", row, (short) monitor.get_Niter()); + bank.setShort("orientation", row, (short) monitor.get_orientation()); + bank.setShort("indicator", row, (short) monitor.get_indicator()); + bank.setShort("status", row, (short) monitor.get_status()); + + row++; + } + } + + return bank; + } + } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java index 522f04f534..30e4fd4772 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java @@ -111,7 +111,7 @@ public boolean init() { this.getConstantsManager().setVariation("default"); - this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction"); + this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction", "AHDC::kftrack:mon"); return true; } @@ -273,6 +273,7 @@ public int compare(Hit a1, Hit a2) { DataBank recoTracksBank = writer.fillAHDCTrackBank(event, AHDC_Tracks); DataBank recoKFTracksBank = writer.fillAHDCKFTrackBank(event, AHDC_Tracks); DataBank AIPredictionBanks = writer.fillAIPrediction(event, predictions); + DataBank recoKFMonitorBank = writer.fillKFMonitorBank(event, AHDC_Tracks); event.appendBank(recoHitsBank); event.appendBank(recoPreClusterBank); @@ -280,6 +281,7 @@ public int compare(Hit a1, Hit a2) { event.appendBank(recoTracksBank); event.appendBank(recoKFTracksBank); event.appendBank(AIPredictionBanks); + event.appendBank(recoKFMonitorBank); if (simulation) { DataBank recoMCBank = writer.fillAHDCMCTrackBank(event); From 174a88153456e300bb56efc05a3be0265a4d73a4 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 12 Nov 2025 04:25:34 -0500 Subject: [PATCH 07/34] optional routine to use mc bank as input of the Kalman filter --- .../jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index cc768f0112..2894ce2b8d 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 @@ -77,18 +77,25 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Initialization material map HashMap materialHashMap = materialGeneration(); int trackId = 0; + DataBank mcBank = event.getBank("MC::Particle"); for (Track track : tracks) { trackId++; track.set_trackId(trackId); // Initialization State Vector - final double x0 = 0.0; - final double y0 = 0.0; - final double z0 = track.get_Z0(); + /*final*/ double x0 = 0.0; + /*final*/ double y0 = 0.0; + /*final*/ double z0 = track.get_Z0(); //final double px0 = track.get_px(); //final double py0 = track.get_py(); - final double pz0 = track.get_pz(); + /*final*/ double pz0 = track.get_pz(); + if (IsMC) { + z0 = mcBank.getFloat("vz", trackId-1)*10; + px0 = mcBank.getFloat("px", trackId-1)*1000; + py0 = mcBank.getFloat("py", trackId-1)*1000; + pz0 = mcBank.getFloat("pz", trackId-1)*1000; + } //final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit From 6452cb1d33894852dafc722ee9b60688ec38f16f Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 12 Nov 2025 09:32:59 -0500 Subject: [PATCH 08/34] save error caovariance in AHDC::kftrack:mon --- etc/bankdefs/hipo4/alert.json | 8 +++++- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 8 ++++++ .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 4 +++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 28 ++++++++++--------- .../org/jlab/rec/ahdc/Track/KFMonitor.java | 8 ++++-- 5 files changed, 40 insertions(+), 16 deletions(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index f7aa7f69b3..dab95e5e4a 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -431,7 +431,13 @@ {"name": "z", "type": "F", "info": "z position of the state (mm)"}, {"name": "px", "type": "F", "info": "px momentum of the state (MeV)"}, {"name": "py", "type": "F", "info": "py momentum of the state (MeV)"}, - {"name": "pz", "type": "F", "info": "pz momentum of the state (MeV)"} + {"name": "pz", "type": "F", "info": "pz momentum of the state (MeV)"}, + {"name": "var_x", "type": "F", "info": "error variance on x (mm^2)"}, + {"name": "var_y", "type": "F", "info": "error variance on y (mm^2)"}, + {"name": "var_z", "type": "F", "info": "error variance on z (mm^2)"}, + {"name": "var_px", "type": "F", "info": "error variance on px (MeV^2)"}, + {"name": "var_py", "type": "F", "info": "error variance on py (MeV^2)"}, + {"name": "var_pz", "type": "F", "info": "error variance on pz (MeV^2)"} ] } ] 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 b6d9c37e6a..a7885abd3d 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 @@ -9,6 +9,7 @@ import org.jlab.rec.ahdc.Track.Track; import org.jlab.rec.ahdc.Track.KFMonitor; import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.linear.RealMatrix; import java.util.ArrayList; @@ -197,6 +198,7 @@ public DataBank fillKFMonitorBank(DataEvent event, ArrayList tracks) { if (track == null) continue; for (KFMonitor monitor : track.get_ListOfKFMonitors()) { RealVector state = monitor.get_state(); + RealMatrix errorCovarianceMatrix = monitor.get_errorCovarianceMatrix(); bank.setFloat("x", row, (float) state.getEntry(0)); bank.setFloat("y", row, (float) state.getEntry(1)); bank.setFloat("z", row, (float) state.getEntry(2)); @@ -208,6 +210,12 @@ public DataBank fillKFMonitorBank(DataEvent event, ArrayList tracks) { bank.setShort("orientation", row, (short) monitor.get_orientation()); bank.setShort("indicator", row, (short) monitor.get_indicator()); bank.setShort("status", row, (short) monitor.get_status()); + bank.setFloat("var_x", row, (float) errorCovarianceMatrix.getEntry(0,0)); + bank.setFloat("var_y", row, (float) errorCovarianceMatrix.getEntry(1,1)); + bank.setFloat("var_z", row, (float) errorCovarianceMatrix.getEntry(2,2)); + bank.setFloat("var_px", row, (float) errorCovarianceMatrix.getEntry(3,3)); + bank.setFloat("var_py", row, (float) errorCovarianceMatrix.getEntry(4,4)); + bank.setFloat("var_pz", row, (float) errorCovarianceMatrix.getEntry(5,5)); row++; } 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 1a184bf5dc..9db6e3c501 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 @@ -320,6 +320,10 @@ public RealVector getStateEstimationVector() { return stateEstimation.copy(); } + public RealMatrix getErrorCovarianceMatrix() { + return errorCovariance.copy(); + } + public double getMomentum() { return Math.sqrt(stateEstimation.getEntry(3) * stateEstimation.getEntry(3) + stateEstimation.getEntry(4) * stateEstimation.getEntry(4) + stateEstimation.getEntry(5) * stateEstimation.getEntry(5)); } 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 2894ce2b8d..a6123c746b 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 @@ -77,7 +77,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Initialization material map HashMap materialHashMap = materialGeneration(); int trackId = 0; - DataBank mcBank = event.getBank("MC::Particle"); + //DataBank mcBank = event.getBank("MC::Particle"); for (Track track : tracks) { trackId++; track.set_trackId(trackId); @@ -90,12 +90,12 @@ private void propagation(ArrayList tracks, DataEvent event, final double //final double py0 = track.get_py(); /*final*/ double pz0 = track.get_pz(); - if (IsMC) { - z0 = mcBank.getFloat("vz", trackId-1)*10; - px0 = mcBank.getFloat("px", trackId-1)*1000; - py0 = mcBank.getFloat("py", trackId-1)*1000; - pz0 = mcBank.getFloat("pz", trackId-1)*1000; - } + //if (IsMC) { + //z0 = mcBank.getFloat("vz", trackId-1)*10; + //px0 = mcBank.getFloat("px", trackId-1)*1000; + //py0 = mcBank.getFloat("py", trackId-1)*1000; + //pz0 = mcBank.getFloat("pz", trackId-1)*1000; + //} //final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit @@ -169,14 +169,16 @@ private void propagation(ArrayList tracks, DataEvent event, final double RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); TrackFitter.setVertexDefined(IsVtxDefined); - + + // save initial state and error covariance matrix + track.add_KFMonitor(new KFMonitor(trackId, -1, -1, -1, -1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); for (int k = 0; k < Niter; k++) { //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: //TrackFitter.ResetErrorCovariance(initialErrorCovariance); for (Indicator indicator : forwardIndicators) { TrackFitter.predict(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { if( k==0 && indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ @@ -184,17 +186,17 @@ private void propagation(ArrayList tracks, DataEvent event, final double } } TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } //System.out.println("--------- BackWard propagation !! ---------"); for (Indicator indicator : backwardIndicators) { TrackFitter.predict(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } } @@ -207,7 +209,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); for (Indicator indicator : forwardIndicators) { PostFitPropagator.predict(indicator); - track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector())); + track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { if( indicator.hit.getHitIdx()>0){ for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java index 190f9cbf1d..f96cb33d4a 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java @@ -1,6 +1,7 @@ package org.jlab.rec.ahdc.Track; import org.apache.commons.math3.linear.RealVector; +import org.apache.commons.math3.linear.RealMatrix; public class KFMonitor { private int trackid; //< trackid @@ -9,15 +10,17 @@ public class KFMonitor { private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]] private int status; //< state just after: prediction (0) or correction (1) private RealVector state; //< vector containing x, y, z, px, py, pz - + private RealMatrix errorCovarianceMatrix; // error covariance matrix, only the diagonal matter (x-x')^2, (y-y')^2, ..., (pz-pz')^2 + // constructor - public KFMonitor(int _trackid, int _Niter, int _orientation, int _indicator, int _status, RealVector _state) { + public KFMonitor(int _trackid, int _Niter, int _orientation, int _indicator, int _status, RealVector _state, RealMatrix _errorCovarianceMatrix) { this.trackid = _trackid; this.Niter = _Niter; this.orientation = _orientation; this.indicator = _indicator; this.status = _status; this.state = _state.copy(); + this.errorCovarianceMatrix = _errorCovarianceMatrix.copy(); } public int get_trackid() {return trackid;} @@ -26,4 +29,5 @@ public KFMonitor(int _trackid, int _Niter, int _orientation, int _indicator, int public int get_indicator() {return indicator;} public int get_status() {return status;} public RealVector get_state() {return state.copy();} + public RealMatrix get_errorCovarianceMatrix() {return errorCovarianceMatrix.copy();} } From 3c0945058d781d5de494c963b618e29d18bf43a4 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 12 Nov 2025 10:07:59 -0500 Subject: [PATCH 09/34] fix kfmonitor id for the initial state --- .../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 a6123c746b..34465fb616 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 @@ -171,7 +171,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double TrackFitter.setVertexDefined(IsVtxDefined); // save initial state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, -1, -1, -1, -1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); for (int k = 0; k < Niter; k++) { //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: From b715c5a5185652c6ed3a72fd0ddd4869a7cea5fd Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 20 Nov 2025 22:49:34 +0100 Subject: [PATCH 10/34] add some comments --- .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 5 +- .../jlab/rec/ahdc/KalmanFilter/Indicator.java | 32 ++++++++--- .../rec/ahdc/KalmanFilter/MaterialMap.java | 6 ++ .../rec/ahdc/KalmanFilter/Propagator.java | 56 +++++++++++-------- .../rec/ahdc/KalmanFilter/RungeKutta4.java | 6 +- .../jlab/rec/ahdc/KalmanFilter/Stepper.java | 30 ++++++---- 6 files changed, 90 insertions(+), 45 deletions(-) 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 a578f7da99..c8582bd237 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java @@ -9,8 +9,11 @@ public class Hit_beam extends Hit { double x,y,z; double r,phi; + // A lot of argument are useless + // FOr clarity, we should remove them + //public Hit_beam(double x, double y , double z) {} public Hit_beam(int superLayer, int layer, int wire, int numWire, double doca, double x, double y , double z) { - super(0, 0, 0, 0, new Line3D(x,y,0,x,y,1), 0); + super(0, 0, 0, 0, new Line3D(x,y,0,x,y,1), 0); // just a line parallel to the beam axis, defined by two end points this.x = x; this.y = y; this.z = z; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java index 09f93054fe..3fe85ccf91 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java @@ -2,13 +2,26 @@ import org.jlab.clas.tracking.kalmanfilter.Material; +// Give the relative positions of the main constituents of the detector +// - beamline +// - target (inner face) +// - target (outer face) +// - ordered list of hits (wires) +// - hit(s) in layer 1 +// - hit(s) in layer 2 +// - +// - hit(s) in layer 8 +// It is used to indicate the Stepper(.java) +// - where to go +// - the medium/material it will cross +// - the step length to be used by RungeKutta4 (RK4) public class Indicator { - public double R; - public double h; - public Hit hit; - public boolean direction; - public Material material; + public double R; // radius position of the detector constituents (cylindrical geometry) + public double h; // step length used by RK4 + public Hit hit; // is null when the stepper moves toward target surfaces + public boolean direction; // true is the forward propagation (beamline to last layer) and false is the backward propagation (last layer to beamline) + public Material material; // medium property public Indicator(double r, double h, Hit hit, boolean direction, Material material) { this.R = r; @@ -17,16 +30,21 @@ public Indicator(double r, double h, Hit hit, boolean direction, Material materi this.direction = direction; this.material = material; } - + + // Not well named + // Tell if we have a measurement or not at this indication (indicator) + // - Each AHDC hit (sense wire) provides a measurement, it is the distance (doca, measured distance to the wire) + // - When the stepper move toward the beamline (backward propagation), we also have a measurement (the distance to the beamline is 0) public boolean haveAHit() { boolean res = this.hit != null; - if (this.R == 0 && !direction) res = true; + if (this.R == 0 && !direction) res = true; // beamline during backward propagation return res; } // This method is not self consistent. We should not have hardcoded values // It is just for testing purpose + // Used by KFMonitor public int getUniqueId() { if (hit != null) { return (hit.getSuperLayer()*10 + hit.getLayer())*100 + hit.getWire(); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java index cfc0ff6495..fcfa011552 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java @@ -5,6 +5,12 @@ import java.util.HashMap; +// In KalmanFilter.java, we have a refinition of this material +// We should merge these codes in the sense that KalmanFilter +// should refer to this class +// +// Also this is not exactly the material for ALERT +// To be modified public class MaterialMap { public static HashMap generateMaterials() { 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 8064c250d7..fe03799437 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 @@ -18,10 +18,11 @@ public class Propagator { public Propagator(RungeKutta4 rungeKutta4) { this.RK4 = rungeKutta4; } - + + // Propagate the stepper toward the next indicator void propagate(Stepper stepper, Indicator indicator) { // ------------------------------------------------------------ - final int maxNbOfStep = 10000; + final int maxNbOfStep = 10000; // do not allow more than 10000 steps (very critical cases) final double R = indicator.R; double dMin = Double.MAX_VALUE; @@ -29,7 +30,8 @@ void propagate(Stepper stepper, Indicator indicator) { // System.out.println("R = " + R); // stepper.print(); - + + // We should use a "while" instead of "for" for (int nbStep = 0; nbStep < maxNbOfStep; nbStep++) { double previous_r = stepper.r(); RK4.doOneStep(stepper); @@ -37,42 +39,45 @@ void propagate(Stepper stepper, Indicator indicator) { // stepper.print(); - if (stepper.direction) { - if (r >= R - 2 - 1.5 * stepper.h) stepper.h = 1e-2; - if (indicator.hit != null) { - if (r >= R - 2) { - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); - if (d < dMin) dMin = d; + if (stepper.direction) { // forward propagation + if (r >= R - 2 - 1.5 * stepper.h) stepper.h = 1e-2; // ? + if (indicator.hit != null) { // the indicator is a hit/wire + if (r >= R - 2) { // only when the stepper is 2 mm closer to the layer of the hit + d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); // distance of the stepper with respect to the wire + if (d < dMin) dMin = d; // the distance d should continue to diminish } - if (r >= R + 2 || d > dMin) { + if (r >= R + 2 || d > dMin) { // r should not be bigger than R + // if the distance d starts to grow again, stop the propagation! + // in that case, the good stepper is the previous one! // System.out.println("dMin = " + dMin); break; } - } else { + } else { // the indicator is the target face (inner or outer) if (r >= R) { - break; + break; // break as soon as r >= R, again, may the good stepper is the previous one } } - } else { - if (r <= R + 2 + 1.5 * stepper.h) stepper.h = 1e-2; - if (R < 5) { - if (stepper.h < 1e-4) stepper.h = 1e-4; - if ((previous_r - r) < 0) { + } else { // backward direction + if (r <= R + 2 + 1.5 * stepper.h) stepper.h = 1e-2; // ? + if (R < 5) { // we are now close the target (outer face is at 3.060 mm) + if (stepper.h < 1e-4) stepper.h = 1e-4; // this step is too short // it increases the computing time + if ((previous_r - r) < 0) { // in the backward propagation, we expect r to decrease; if not stop the propagation! break; } } - if (indicator.hit != null) { - if (r <= R + 2) { + if (indicator.hit != null) { // wire/hit or beamline + if (r <= R + 2) { // here R < r and r decreases, we enter here when r becomes 2 mm closer to R + // we compute the distance, we expect it to decrease d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); if (d < dMin) dMin = d; } - if (r <= R - 2 || d > dMin) { + if (r <= R - 2 || d > dMin) { // if the distance grows again, stop the propagation // System.out.println("dMin = " + dMin); break; } - } else { - if (r <= R) { - break; + } else { // the indicator is the target face (inner or outer) + if (r <= R) { + break; // break as soon as r <= R, again, may the good stepper is the previous one } } @@ -88,7 +93,10 @@ public RealVector f(Stepper stepper, Indicator indicator) { propagate(stepper, indicator); return new ArrayRealVector(stepper.y); } - + + // ------------------------------- + // not used, to be deleted + // ------------------------------- public void propagateAndWrite(Stepper stepper, Indicator indicator, Writer writer) { // ------------------------------------------------------------ final int maxNbOfStep = 10000; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java index 8830b959e0..e66f1a3a56 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java @@ -5,6 +5,9 @@ import java.util.Arrays; +// Compute the next trajectory point using the equation of motion +// the path length is given by the stepper +// doneOneStep() is the main method, it is used in the Propagator.propagate() method public class RungeKutta4 { private final int numberOfVariables; private final PDGParticle particle; @@ -96,7 +99,8 @@ public void doOneStep(Stepper stepper) { } private double[] f(double[] y) { - double charge = 1.0; + double charge = 1.0; // the charge should by given by the particle object + // for now, we assume it is a proton double pModuleInverse = 1.0 / Math.sqrt(y[3] * y[3] + y[4] * y[4] + y[5] * y[5]); double k = charge * PhysicsConstants.speedOfLight() * 10 * pModuleInverse; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java index cee2861f9c..1d6b52709c 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java @@ -4,21 +4,25 @@ import java.util.Arrays; +// A trajectory point computed by the Kalman Filter public class Stepper { - public double[] y; - public double h; - public Material material; - public double s; - public double sTot = 0; - public double s_drift = 0; - public boolean is_in_drift = false; - public double dEdx; - public boolean direction; + public double[] y; // position and momentum coordinates + public double h; // + public Material material; + public double s; // path length during the propagation between the current indicator and the next one given in initialize() + public double sTot = 0; // total lenght with respect to the beamline // increase in forward propagation, reach the maximum in the last layer and decrease in the backward propagation // cf doOneStep() method in RungeKutta4.java + public double s_drift = 0; // cf doOneStep() method in RungeKutta4.java + public boolean is_in_drift = false; // cf doOneStep() method in RungeKutta4.java + public double dEdx; // deposited energy during the propagation between the current indicator and the next one given in initialize() + public boolean direction; // true is the forward propagation (beamline to last layer) and false is the backward propagation (last layer to beamline) public Stepper(double[] y) { this.y = Arrays.copyOf(y, y.length); } - + + // The initialisation tells the stepper where to go (i.e toward the `indicator`) + // Cf. notes in Indicator.java + // The initilisation is done before any propagation between indicators, in the predict() method of KFitter.java public void initialize(Indicator indicator) { this.s = 0; this.dEdx = 0; @@ -28,11 +32,13 @@ public void initialize(Indicator indicator) { this.is_in_drift = false; } - + + // radius, distance to the beamline public double r() { return Math.hypot(y[0], y[1]); } - + + // momentum public double p() { return Math.sqrt(y[3] * y[3] + y[4] * y[4] + y[5] * y[5]); } From 8a5843fea554c33861d15a2887d32744332087c0 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 21 Nov 2025 13:27:19 +0100 Subject: [PATCH 11/34] remove everything related to sign ambiguity --- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 37 ----- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 80 ++-------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 146 ++++++------------ 3 files changed, 60 insertions(+), 203 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java index 0aeabb8af0..8e649d3139 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java @@ -23,7 +23,6 @@ public class Hit implements Comparable { private final Line3D line3D_plus; private final Line3D line3D_minus; private int hitidx; - private int hitsign; // Comparison with: common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java // here, SuperLayer, Layer, Wire, start from 1 @@ -37,7 +36,6 @@ public Hit(int superLayer, int layer, int wire, int numWire, Line3D line, double this.numWires = numWire; this.adc = 0;//placeholder this.hitidx = -1; - this.hitsign = 0; this.phi = line.midpoint().vectorFrom(0,0,0).phi(); //System.out.println(" superlayer " + this.superLayer + " layer " + this.layer + " wire " + this.wire + " wx " + wx + " wy " + wy + " wx_end " + wx_end + " wy_end " + wy_end + " phi " + this.phi); @@ -61,27 +59,10 @@ public RealVector get_Vector() { return new ArrayRealVector(new double[]{this.doca}); } - //hit measurement vector in 1 dimension with sign: if sign = 0, return doca, otherwise return 0 - public RealVector get_Vector(int sign, boolean goodsign) { - if(sign == 0 || goodsign){ - return new ArrayRealVector(new double[]{this.doca}); - }else{ - return new ArrayRealVector(new double[]{0.0}); - } - } - public RealMatrix get_MeasurementNoise() { return new Array2DRowRealMatrix(new double[][]{{0.0225}}); } - public RealMatrix get_MeasurementNoise(boolean goodsign) { - if(goodsign){ - return new Array2DRowRealMatrix(new double[][]{{0.0225}}); - }else{ - return new Array2DRowRealMatrix(new double[][]{{2*this.doca*this.doca}}); - } - } - public double doca() { return doca; } @@ -100,16 +81,6 @@ public double distance(Point3D point3D) { return this.line3D.distance(point3D).length(); } - public double distance(Point3D point3D, int sign, boolean goodsign) { - //if(sign!=0) - //System.out.println(" r " + this.r + " phi " + this.phi + " doca " + this.doca + " sign " + sign + " distance " + this.line3D.distance(point3D).length() + " (sign 0) " + this.line3D_plus.distance(point3D).length() + " (sign+) " + this.line3D_minus.distance(point3D).length() + " (sign-) "); - if(!goodsign){ - if(sign>0)return this.line3D_plus.distance(point3D).length(); - if(sign<0)return this.line3D_minus.distance(point3D).length(); - } - return this.line3D.distance(point3D).length(); - } - @Override public int compareTo(Hit o) { System.out.println("r = " + r + " other r = " + o.r); @@ -169,13 +140,5 @@ public void setHitIdx(int idx) { this.hitidx = idx; } - public int getSign() { - return hitsign; - } - - public void setSign(int sign) { - this.hitsign = sign; - } - } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 9db6e3c501..5bc5b1b4de 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 @@ -10,11 +10,11 @@ public class KFitter { private RealVector stateEstimation; - private RealMatrix errorCovariance; + private RealMatrix errorCovariance; public final Stepper stepper; private final Propagator propagator; - public double chi2 = 0; - // masses/energies in MeV + public double chi2 = 0; + // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; private boolean isvertexdefined = false; @@ -73,12 +73,12 @@ public void predict(Indicator indicator) throws Exception { } public void correct(Indicator indicator) { - RealVector z, z_plus, z_minus; + RealVector z, z_plus, z_minus; RealMatrix measurementNoise; RealMatrix measurementMatrix; RealVector h; if (indicator.R == 0.0 && !indicator.direction) { - double z_beam_res_sq = 1.e10;//in mm + double z_beam_res_sq = 1.e10;//in mm if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 2. mm resolution measurementNoise = new Array2DRowRealMatrix( @@ -91,21 +91,10 @@ public void correct(Indicator indicator) { h = h_beam(stateEstimation);//3x1 z = indicator.hit.get_Vector_beam();//0! } else { - //System.out.println(" hit r " + indicator.hit.r() + " hit phi " + indicator.hit.phi() + " phi wire (-zl/2) " + indicator.hit.phi(-150.0) + " phi wire (0) " + indicator.hit.phi(0.0) + " phi wire (+zl/2) " + indicator.hit.phi(150.) + " state x " + stateEstimation.getEntry(0) + " state y " + stateEstimation.getEntry(1) + " state z " + stateEstimation.getEntry(2) ); - boolean goodsign = true; - if(indicator.hit.getSign()!=0){ - double dphi = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0))-indicator.hit.phi(stateEstimation.getEntry(2)); - if(dphi*indicator.hit.getSign()<0)goodsign = false; - //System.out.println(" hit r " + indicator.hit.r() + " phi wire (z) " + indicator.hit.phi(stateEstimation.getEntry(2)) + " phi state " + Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)) + " sign " + indicator.hit.getSign() + " good? " + goodsign ); - } - //measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 - measurementNoise = indicator.hit.get_MeasurementNoise(goodsign);//1x1 - measurementMatrix = H(stateEstimation, indicator);//6x1 - //measurementMatrix = H(stateEstimation, indicator, goodsign);//6x1 - h = h(stateEstimation, indicator);//1x1 - //h = h(stateEstimation, indicator, goodsign);//1x1 + measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 + measurementMatrix = H(stateEstimation, indicator);//6x1 + h = h(stateEstimation, indicator);//1x1 z = indicator.hit.get_Vector();//1x1 - //z = indicator.hit.get_Vector(indicator.hit.getSign(), goodsign);//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); @@ -123,6 +112,7 @@ public void correct(Indicator indicator) { // update estimate with measurement z(k) xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error P(k) = (I - K * H) * P(k)- + // The Joseph's form is numerically more stable P(k) = (I - K * H) * P(k)- (I - K * H)' - H*R*K' RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); // Numerically more stable !! RealMatrix tmpMatrix = identity.subtract(kalmanGain.multiply(measurementMatrix)); @@ -137,20 +127,9 @@ public double residual(Indicator indicator) { return indicator.hit.doca()-d; } - //function for left-right disambiguation - public int wire_sign(Indicator indicator) {//let's decide: positive when (phi state - phi wire) > 0 - double phi_state = Math.atan2(stateEstimation.getEntry(1), stateEstimation.getEntry(0)); - double phi_wire = indicator.hit.phi(stateEstimation.getEntry(2)); - if( (phi_state-phi_wire)/Math.abs(phi_state-phi_wire)>0 ){ - return +1; - }else{ - return -1; - } - } - - public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ - this.errorCovariance = initialErrorCovariance; - } + public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ + this.errorCovariance = initialErrorCovariance; + } private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { @@ -192,12 +171,6 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { //measurement matrix in 1 dimension: minimize distance - doca private RealVector h(RealVector x, Indicator indicator) { double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - //double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign()); - return MatrixUtils.createRealVector(new double[]{d}); - } - - private RealVector h(RealVector x, Indicator indicator, boolean goodsign) { - double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2)), indicator.hit.getSign(), goodsign); return MatrixUtils.createRealVector(new double[]{d}); } @@ -230,35 +203,6 @@ private RealMatrix H(RealVector x, Indicator indicator) { return (doca_plus - doca_minus) / (2 * h); } - //measurement matrix in 1 dimension: minimize distance - doca - private RealMatrix H(RealVector x, Indicator indicator, boolean goodsign) { - - double ddocadx = subfunctionH(x, indicator, 0, goodsign); - double ddocady = subfunctionH(x, indicator, 1, goodsign); - double ddocadz = subfunctionH(x, indicator, 2, goodsign); - double ddocadpx = subfunctionH(x, indicator, 3, goodsign); - double ddocadpy = subfunctionH(x, indicator, 4, goodsign); - double ddocadpz = subfunctionH(x, indicator, 5, goodsign); - - // As per my understanding: ddocadx,y,z -> = dr/dx,y,z, etc - return MatrixUtils.createRealMatrix(new double[][]{ - {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); - } - - double subfunctionH(RealVector x, Indicator indicator, int i, boolean goodsign) { - double h = 1e-8;// in mm - RealVector x_plus = x.copy(); - RealVector x_minus = x.copy(); - - x_plus.setEntry(i, x_plus.getEntry(i) + h); - x_minus.setEntry(i, x_minus.getEntry(i) - h); - - double doca_plus = h(x_plus, indicator, goodsign).getEntry(0); - double doca_minus = h(x_minus, indicator, goodsign).getEntry(0); - - return (doca_plus - doca_minus) / (2 * h); - } - private RealVector h_beam(RealVector x) { double xx = x.getEntry(0); 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 34465fb616..306b60d6d2 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java @@ -34,125 +34,71 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 10; - private final boolean IsVtxDefined = false; + private final int Niter = 10; // number of iterations of the Kalman Filter + private final boolean IsVtxDefined = false; // never used so far + // can be useful to take into account the vertex of the electron private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { try { - double vz_constraint = 0; - if(IsMC) {//If simulation read MC::Particle Bank ------------------------------------------------ - DataBank bankParticle = event.getBank("MC::Particle"); - double vxmc = bankParticle.getFloat("vx", 0)*10;//mm - double vymc = bankParticle.getFloat("vy", 0)*10;//mm - double vzmc = bankParticle.getFloat("vz", 0)*10;//mm - double pxmc = bankParticle.getFloat("px", 0)*1000;//MeV - double pymc = bankParticle.getFloat("py", 0)*1000;//MeV - double pzmc = bankParticle.getFloat("pz", 0)*1000;//MeV - double p_mc = java.lang.Math.sqrt(pxmc*pxmc+pymc*pymc+pzmc*pzmc); - //System.out.println("MC track: vz: " + vzmc*10 + " px: " + pxmc*1000 + " py: " + pymc*1000 + " pz: " + pzmc*1000 + "; p = " + p_mc*1000);//convert p to MeV, v to mm - - ArrayList sim_hits = new ArrayList<>(); - sim_hits.add(new Point3D(0, 0, vzmc)); - - DataBank bankMC = event.getBank("MC::True"); - for (int i = 0; i < bankMC.rows(); i++) { - if (bankMC.getInt("pid", i) == 2212) { - float x = bankMC.getFloat("avgX", i); - float y = bankMC.getFloat("avgY", i); - float z = bankMC.getFloat("avgZ", i); - // System.out.println("r_sim = " + Math.hypot(x, y)); - sim_hits.add(new Point3D(x, y, z)); - } - } - vz_constraint = vzmc; - } + double vz_constraint = 0; // to be linked to the electron vertex // Initialization --------------------------------------------------------------------- final PDGParticle proton = PDGDatabase.getParticleById(2212); final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; + + //DataBank mcBank = event.getBank("MC::Particle"); // Initialization material map HashMap materialHashMap = materialGeneration(); + // Loop over tracks + // Each track candidate form the HelixFitter is independently processed by the Kalman Filter int trackId = 0; - //DataBank mcBank = event.getBank("MC::Particle"); for (Track track : tracks) { trackId++; track.set_trackId(trackId); // Initialization State Vector - /*final*/ double x0 = 0.0; - /*final*/ double y0 = 0.0; - /*final*/ double z0 = track.get_Z0(); - //final + double x0 = 0.0; + double y0 = 0.0; + double z0 = track.get_Z0(); double px0 = track.get_px(); - //final double py0 = track.get_py(); - /*final*/ double pz0 = track.get_pz(); - //if (IsMC) { - //z0 = mcBank.getFloat("vz", trackId-1)*10; - //px0 = mcBank.getFloat("px", trackId-1)*1000; - //py0 = mcBank.getFloat("py", trackId-1)*1000; - //pz0 = mcBank.getFloat("pz", trackId-1)*1000; - //} - //final double p_init = java.lang.Math.sqrt(px0*px0+py0*py0+pz0*pz0); + double pz0 = track.get_pz(); + // using or not mc + // will be deleted in the final code + /*if (IsMC) { + z0 = mcBank.getFloat("vz", trackId-1)*10; + px0 = mcBank.getFloat("px", trackId-1)*1000; + py0 = mcBank.getFloat("py", trackId-1)*1000; + pz0 = mcBank.getFloat("pz", trackId-1)*1000; + }*/ double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit + // Here the hits that will be used by the Kalman Filter are initialised from the one AHDC::hits bank + // Why defining a new type of hit ? + // Unifying could be useful for the processing time. To be checked later. ArrayList AHDC_hits = track.getHits(); ArrayList KF_hits = new ArrayList<>(); double ADCTot = 0; - //System.out.println(" px " + y[3] + " py " + y[4] +" pz " + y[5] +" vz " + y[2] + " number of hits: " + AHDC_hits.size() + " MC hits? " + sim_hits.size()); track.set_n_hits(AHDC_hits.size()); for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { - Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getLine(), AHDC_hit.getDoca()); - hit.setADC(AHDC_hit.getADC()); - hit.setHitIdx(AHDC_hit.getId()); - hit.setSign(0); - ADCTot+=AHDC_hit.getADC(); - // set track id - AHDC_hit.setTrackId(trackId); - - //System.out.println( " r = " + hit.r() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca() ); - // Do delete hit with same radius - boolean phi_rollover = false; - boolean aleardyHaveR = false; - for (Hit o: KF_hits){ - if (o.r() == hit.r()){ - aleardyHaveR = true; - // //sign+ means (phi track - phi wire) > 0 - // if(o.phi()>hit.phi()){ - // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - // o.setSign(-1); - // hit.setSign(+1); - // }else{ - // phi_rollover = true; - // hit.setSign(-1); - // o.setSign(+1); - // } - // }else{ - // if(Math.abs(o.phi()-hit.phi())< 2*Math.toRadians(360./o.getNumWires()) ){ - // hit.setSign(-1); - // o.setSign(+1); - // }else{ - // phi_rollover = true; - // o.setSign(-1); - // hit.setSign(+1); - // } - // } - // //System.out.println( " r = " + o.r() + " o.phi = " + o.phi() + " o.doca = " + o.getDoca()*o.getSign() + " hit.phi " + hit.phi() +" hit.doca = " + hit.getDoca()*hit.getSign() + " angle between wires: " + Math.toRadians(360./hit.getNumWires()) + " >= ? angle covered by docas: " + Math.atan( (o.getDoca()+hit.getDoca())/o.r() ) ); - } - } - if(!aleardyHaveR)KF_hits.add(hit); - // if (phi_rollover){ - // KF_hits.add(KF_hits.size()-1, hit); - // }else{ - // KF_hits.add(hit); - // } + // Define ahdc.KalmanFilter.Hit from ahdc.Hit.Hit + Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getLine(), AHDC_hit.getDoca()); + hit.setADC(AHDC_hit.getADC()); + hit.setHitIdx(AHDC_hit.getId()); + ADCTot+=AHDC_hit.getADC(); + // set track id + AHDC_hit.setTrackId(trackId); + KF_hits.add(hit); } + // add a routine to sort the list of hits double zbeam = 0; if(IsVtxDefined)zbeam = vz_constraint;//test + // Define forward and backward indicator + // cf. Indicator.java final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam); @@ -165,28 +111,32 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Initialization of the Kalman Fitter RealVector initialStateEstimate = new ArrayRealVector(stepper.y); - //first 3 lines in cm^2; last 3 lines in MeV^2 + //first 3 lines in mm^2; last 3 lines in MeV^2 RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); TrackFitter.setVertexDefined(IsVtxDefined); - // save initial state and error covariance matrix + // KFmonitor: save initial state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); for (int k = 0; k < Niter; k++) { //System.out.println("--------- ForWard propagation !! ---------"); //Reset error covariance: - //TrackFitter.ResetErrorCovariance(initialErrorCovariance); + //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed for (Indicator indicator : forwardIndicators) { + // Prediction TrackFitter.predict(indicator); + // KFMonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { - if( k==0 && indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(TrackFitter.residual(indicator)); - } - } - TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + // I don't see the utility of this + if( k==0 && indicator.hit.getHitIdx()>0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(TrackFitter.residual(indicator)); + } + } + // Correction only if we have a measure (hit) + TrackFitter.correct(indicator); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } From e8c1fa56ea9ca21554816ecda2ca89f7837d1f8c Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 21 Nov 2025 13:34:36 +0100 Subject: [PATCH 12/34] fix indentation --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 82 +++++++++---------- 1 file changed, 40 insertions(+), 42 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 306b60d6d2..8fa4c78f4e 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 @@ * TODO : - Fix multi hit on the same layer * - Optimize measurement noise and probably use doca as weight * - Fix the wire number (-1) - * - Iterate thought multiple tracks per event * - use a flag for simulation * - add target in the map * - move map to initialization engine @@ -118,37 +117,38 @@ private void propagation(ArrayList tracks, DataEvent event, final double // KFmonitor: save initial state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + // Loop over number of iterations for (int k = 0; k < Niter; k++) { - //System.out.println("--------- ForWard propagation !! ---------"); - //Reset error covariance: - //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed - for (Indicator indicator : forwardIndicators) { - // Prediction - TrackFitter.predict(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - // I don't see the utility of this - if( k==0 && indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(TrackFitter.residual(indicator)); + //System.out.println("--------- ForWard propagation !! ---------"); + //Reset error covariance: + //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed + for (Indicator indicator : forwardIndicators) { + // Prediction + TrackFitter.predict(indicator); + // KFMonitor: save state and error covariance matrix + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + if (indicator.haveAHit()) { + // I don't see the utility of this + if( k==0 && indicator.hit.getHitIdx()>0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(TrackFitter.residual(indicator)); + } } + // Correction only if we have a measure (hit) + TrackFitter.correct(indicator); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } - // Correction only if we have a measure (hit) - TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - } - } + } - //System.out.println("--------- BackWard propagation !! ---------"); - for (Indicator indicator : backwardIndicators) { - TrackFitter.predict(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - TrackFitter.correct(indicator); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - } - } + //System.out.println("--------- BackWard propagation !! ---------"); + for (Indicator indicator : backwardIndicators) { + TrackFitter.predict(indicator); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + if (indicator.haveAHit()) { + TrackFitter.correct(indicator); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + } + } } @@ -158,15 +158,15 @@ private void propagation(ArrayList tracks, DataEvent event, final double //Residual, path and AHDC exit momentum calculation post fit: KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); for (Indicator indicator : forwardIndicators) { - PostFitPropagator.predict(indicator); - track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - if( indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(PostFitPropagator.residual(indicator)); - } - } - } + PostFitPropagator.predict(indicator); + track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + if (indicator.haveAHit()) { + if( indicator.hit.getHitIdx()>0){ + for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ + if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(PostFitPropagator.residual(indicator)); + } + } + } } double s = PostFitPropagator.stepper.sTot; @@ -179,9 +179,9 @@ private void propagation(ArrayList tracks, DataEvent event, final double double sum_residuals = 0; double chi2 = 0; for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { - sum_adc += AHDC_hit.getADC(); - sum_residuals += AHDC_hit.getResidual(); - chi2 += Math.pow(AHDC_hit.getResidual(),2.0); + sum_adc += AHDC_hit.getADC(); + sum_residuals += AHDC_hit.getResidual(); + chi2 += Math.pow(AHDC_hit.getResidual(),2.0); } track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); @@ -193,8 +193,6 @@ private void propagation(ArrayList tracks, DataEvent event, final double } catch (Exception e) { // e.printStackTrace(); } - - } private HashMap materialGeneration() { From 6b459f378194ebba0f9dbc8cf46b0ae7bbdb87b9 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 21 Nov 2025 15:40:35 +0100 Subject: [PATCH 13/34] unify ahdc.Hit.Hit and ahdc.KalmanFilter.Hit --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 23 ++- .../org/jlab/rec/ahdc/KalmanFilter/Hit.java | 144 ------------------ .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 18 ++- .../jlab/rec/ahdc/KalmanFilter/Indicator.java | 7 +- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 10 +- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 68 ++++----- 6 files changed, 69 insertions(+), 201 deletions(-) delete mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.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 26007ae847..f419517fe3 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 @@ -4,6 +4,11 @@ import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; +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; + public class Hit implements Comparable { private final double thster = Math.toRadians(20.0); @@ -45,7 +50,7 @@ public void setWirePosition(AlertDCDetector factory) { //System.out.println(" superlayer " + this.superLayerId + " layer " + this.layerId + " wire " + this.wireId + " R_layer " + R_layer + " wx " + wx + " wy " + wy); wireLine = factory.getSector(1).getSuperlayer(superLayerId).getLayer(layerId).getComponent(wireId).getLine(); Point3D end = wireLine.end(); - this.nbOfWires = factory.getSector(1).getSuperlayer(superLayerId).getLayer(layerId).getNumComponents(); + this.nbOfWires = factory.getSector(1).getSuperlayer(superLayerId).getLayer(layerId).getNumComponents(); this.phi = end.vectorFrom(0, 0, 0).phi(); this.radius = end.distance(0, 0, end.z()); this.x = end.x(); @@ -146,4 +151,20 @@ public void setTrackId(int _trackId) { this.trackId = _trackId; } + public RealVector get_Vector() { + return new ArrayRealVector(new double[]{this.doca}); + } + + public RealMatrix get_MeasurementNoise() { + return new Array2DRowRealMatrix(new double[][]{{0.0225}}); + } + + 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.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java deleted file mode 100644 index 8e649d3139..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit.java +++ /dev/null @@ -1,144 +0,0 @@ -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.Plane3D; -import org.jlab.geom.prim.Point3D; -import org.jlab.geom.prim.Vector3D; - -public class Hit implements Comparable { - - private final int superLayer; - private final int layer; - private final int wire; - private final double r; - private final double phi; - private final double doca; - private double adc; - private final double numWires; - private final Line3D line3D; - private final Line3D line3D_plus; - private final Line3D line3D_minus; - private int hitidx; - - // Comparison with: common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java - // here, SuperLayer, Layer, Wire, start from 1 - // in AlertDCFactory, same variables start from 1 - public Hit(int superLayer, int layer, int wire, int numWire, Line3D line, double doca) { - this.superLayer = superLayer; - this.layer = layer; - this.wire = wire; - this.r = line.end().distance(0, 0, line.end().z()); - this.doca = doca; - this.numWires = numWire; - this.adc = 0;//placeholder - this.hitidx = -1; - - this.phi = line.midpoint().vectorFrom(0,0,0).phi(); - //System.out.println(" superlayer " + this.superLayer + " layer " + this.layer + " wire " + this.wire + " wx " + wx + " wy " + wy + " wx_end " + wx_end + " wy_end " + wy_end + " phi " + this.phi); - - this.line3D = line; - - //calculate the "virtual" left and right wires accounting for the DOCA - double deltaphi = Math.asin(this.doca/r); - this.line3D_plus = new Line3D(); - this.line3D_plus.copy(line); - this.line3D_plus.rotateZ(deltaphi); - - this.line3D_minus = new Line3D(); - this.line3D_minus.copy(line); - this.line3D_minus.rotateZ(deltaphi); - - } - - //hit measurement vector in 1 dimension: minimize distance - doca - public RealVector get_Vector() { - return new ArrayRealVector(new double[]{this.doca}); - } - - public RealMatrix get_MeasurementNoise() { - return new Array2DRowRealMatrix(new double[][]{{0.0225}}); - } - - public double doca() { - return doca; - } - - public double r() {return r;} - - public double phi() {return phi;}//at z = 0; - - public double phi(double z) { - return this.line3D.lerpPoint((z-line3D.origin().z())/line3D.length()).vectorFrom(0,0,0).phi(); - } - - public Line3D line() {return line3D;} - - public double distance(Point3D point3D) { - return this.line3D.distance(point3D).length(); - } - - @Override - public int compareTo(Hit o) { - System.out.println("r = " + r + " other r = " + o.r); - return Double.compare(r, o.r); - } - - @Override - public String toString() { - return "Hit{" + "superLayer=" + superLayer + ", layer=" + layer + ", wire=" + wire + ", r=" + r + ", doca=" + doca + '}'; - } - - public RealVector get_Vector_beam() { - return null; - } - - public int getSuperLayer() { - return superLayer; - } - - public int getLayer() { - return layer; - } - - public int getWire() { - return wire; - } - - public double getR() { - return r; - } - - public double getDoca() { - return doca; - } - - public double getADC() { - return adc; - } - - public void setADC(double _adc) { - this.adc = _adc; - } - - public Line3D getLine3D() { - return line3D; - } - - public double getNumWires() { - return numWires; - } - - public int getHitIdx() { - return hitidx; - } - - public void setHitIdx(int idx) { - this.hitidx = idx; - } - -} - 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 c8582bd237..4e5fc3be09 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 @@ -3,25 +3,35 @@ import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Line3D; +import org.jlab.rec.ahdc.Hit.Hit; +// A weird object that just want to be considered as a Hit +// but it is not! +// It is due to the notion of Indicator that want to unify wire hits, beamline or target surfaces public class Hit_beam extends Hit { - double x,y,z; - double r,phi; + private double x,y,z; + private double r,phi; + Line3D beamline; // A lot of argument are useless // FOr clarity, we should remove them //public Hit_beam(double x, double y , double z) {} - public Hit_beam(int superLayer, int layer, int wire, int numWire, double doca, double x, double y , double z) { - super(0, 0, 0, 0, new Line3D(x,y,0,x,y,1), 0); // just a line parallel to the beam axis, defined by two end points + public Hit_beam(double x, double y , double z) { + super(-1, -1, -1, -1, 0, 0, -1); // just a line parallel to the beam axis, defined by two end points this.x = x; this.y = y; this.z = z; this.r = Math.hypot(x,y); this.phi = Math.atan2(y,x); + beamline = new Line3D(x,y,0,x,y,1); } public RealVector get_Vector_beam() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } + + public Line3D getLine() { + return beamline; + } } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java index 3fe85ccf91..d4f88e55d7 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java @@ -1,6 +1,7 @@ package org.jlab.rec.ahdc.KalmanFilter; import org.jlab.clas.tracking.kalmanfilter.Material; +import org.jlab.rec.ahdc.Hit.Hit; // Give the relative positions of the main constituents of the detector // - beamline @@ -43,11 +44,11 @@ public boolean haveAHit() { // This method is not self consistent. We should not have hardcoded values - // It is just for testing purpose - // Used by KFMonitor + // It is just for testing + // Only used by KFMonitor public int getUniqueId() { if (hit != null) { - return (hit.getSuperLayer()*10 + hit.getLayer())*100 + hit.getWire(); + return (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(); } else if (Math.abs(R) < 1e-9) { // beamline 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 5bc5b1b4de..aa682d4fea 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java @@ -73,7 +73,7 @@ public void predict(Indicator indicator) throws Exception { } public void correct(Indicator indicator) { - RealVector z, z_plus, z_minus; + RealVector z; RealMatrix measurementNoise; RealMatrix measurementMatrix; RealVector h; @@ -124,7 +124,7 @@ public void correct(Indicator indicator) { public double residual(Indicator indicator) { double d = indicator.hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); - return indicator.hit.doca()-d; + return indicator.hit.getDoca()-d; } public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ @@ -249,12 +249,6 @@ private RealMatrix H_beam(RealVector x) { }); } - - private RealVector innovation(RealVector z, Indicator indicator) { - - return z.subtract(h(stateEstimation, indicator)); - } - /** * Returns a copy of the current state estimation vector. * 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 8fa4c78f4e..e0dbb8f6dc 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 @@ -12,6 +12,7 @@ import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Track.Track; +import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.KFMonitor; import java.util.ArrayList; @@ -75,31 +76,17 @@ private void propagation(ArrayList tracks, DataEvent event, final double }*/ double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit - // Here the hits that will be used by the Kalman Filter are initialised from the one AHDC::hits bank - // Why defining a new type of hit ? - // Unifying could be useful for the processing time. To be checked later. - ArrayList AHDC_hits = track.getHits(); - ArrayList KF_hits = new ArrayList<>(); - double ADCTot = 0; - track.set_n_hits(AHDC_hits.size()); - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { - // Define ahdc.KalmanFilter.Hit from ahdc.Hit.Hit - Hit hit = new Hit(AHDC_hit.getSuperLayerId(), AHDC_hit.getLayerId(), AHDC_hit.getWireId(), AHDC_hit.getNbOfWires(), AHDC_hit.getLine(), AHDC_hit.getDoca()); - hit.setADC(AHDC_hit.getADC()); - hit.setHitIdx(AHDC_hit.getId()); - ADCTot+=AHDC_hit.getADC(); - // set track id - AHDC_hit.setTrackId(trackId); - KF_hits.add(hit); - } + ArrayList AHDC_hits = track.getHits(); + //----------------------------------------- // add a routine to sort the list of hits + //----------------------------------------- double zbeam = 0; if(IsVtxDefined)zbeam = vz_constraint;//test // Define forward and backward indicator // cf. Indicator.java - final ArrayList forwardIndicators = forwardIndicators(KF_hits, materialHashMap); - final ArrayList backwardIndicators = backwardIndicators(KF_hits, materialHashMap, zbeam); + final ArrayList forwardIndicators = forwardIndicators(AHDC_hits, materialHashMap); + final ArrayList backwardIndicators = backwardIndicators(AHDC_hits, materialHashMap, zbeam); // Start propagation Stepper stepper = new Stepper(y); @@ -129,13 +116,12 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { // I don't see the utility of this - if( k==0 && indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidualPrefit(TrackFitter.residual(indicator)); - } + if( k==0 && indicator.hit.getId()>0){ // first iteration and indicator != beamline (because its id is -1) + indicator.hit.setResidualPrefit(TrackFitter.residual(indicator)); } // Correction only if we have a measure (hit) TrackFitter.correct(indicator); + // KFMonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } @@ -143,9 +129,11 @@ private void propagation(ArrayList tracks, DataEvent event, final double //System.out.println("--------- BackWard propagation !! ---------"); for (Indicator indicator : backwardIndicators) { TrackFitter.predict(indicator); + // KFMonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { TrackFitter.correct(indicator); + // KFMonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } @@ -159,36 +147,34 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); for (Indicator indicator : forwardIndicators) { PostFitPropagator.predict(indicator); + // KFMonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { - if( indicator.hit.getHitIdx()>0){ - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits){ - if(AHDC_hit.getId()==indicator.hit.getHitIdx())AHDC_hit.setResidual(PostFitPropagator.residual(indicator)); - } + if( indicator.hit.getId()>0){ + indicator.hit.setResidual(TrackFitter.residual(indicator)); } } } + // Fill track and hit bank double s = PostFitPropagator.stepper.sTot; - double dEdx = ADCTot / s; double p_drift = PostFitPropagator.stepper.p(); - - // At this stage, all relevants AHDC hits are filled - // Compute sum_adc, sum_residuals and chi2 int sum_adc = 0; double sum_residuals = 0; double chi2 = 0; - for (org.jlab.rec.ahdc.Hit.Hit AHDC_hit : AHDC_hits) { - sum_adc += AHDC_hit.getADC(); - sum_residuals += AHDC_hit.getResidual(); - chi2 += Math.pow(AHDC_hit.getResidual(),2.0); + for (Hit hit : AHDC_hits) { + hit.setTrackId(trackId); + sum_adc += hit.getADC(); + sum_residuals += hit.getResidual(); + chi2 += Math.pow(hit.getResidual(),2.0); } track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); track.set_chi2(chi2); track.set_p_drift_kf(p_drift); - track.set_dEdx_kf(dEdx); + track.set_dEdx_kf(sum_adc/s); track.set_path_kf(s); + track.set_n_hits(AHDC_hits.size()); }//end of loop on track candidates } catch (Exception e) { // e.printStackTrace(); @@ -250,7 +236,7 @@ ArrayList forwardIndicators(ArrayList hitArrayList, HashMap backwardIndicators(ArrayList hitArrayList, HashMap backwardIndicators = new ArrayList<>(); //R, h, defined in mm! for (int i = hitArrayList.size() - 2; i >= 0; i--) { - backwardIndicators.add(new Indicator(hitArrayList.get(i).r(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); + backwardIndicators.add(new Indicator(hitArrayList.get(i).getRadius(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); } backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); - Hit hit = new Hit_beam(0, 0, 0, 0, 0, 0, 0, 0); + Hit hit = new Hit_beam(0, 0, 0); backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); return backwardIndicators; } @@ -272,11 +258,11 @@ ArrayList backwardIndicators(ArrayList hitArrayList, HashMap backwardIndicators = new ArrayList<>(); //R, h, defined in mm! for (int i = hitArrayList.size() - 2; i >= 0; i--) { - backwardIndicators.add(new Indicator(hitArrayList.get(i).r(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); + backwardIndicators.add(new Indicator(hitArrayList.get(i).getRadius(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); } backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); - Hit hit = new Hit_beam(0, 0, 0, 0, 0, 0, 0, vz); + Hit hit = new Hit_beam(0, 0, vz); backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); return backwardIndicators; } From 58ad6c583dbc260704475c5837dcf474277353c3 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Fri, 21 Nov 2025 17:11:18 +0100 Subject: [PATCH 14/34] fix indentation --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 19 +++++++++---------- .../java/org/jlab/rec/ahdc/Hit/HitReader.java | 6 +++--- 2 files changed, 12 insertions(+), 13 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 f419517fe3..e9e9d84115 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 @@ -11,7 +11,6 @@ public class Hit implements Comparable { - private final double thster = Math.toRadians(20.0); private final int id; private final int superLayerId; private final int layerId; @@ -20,7 +19,7 @@ public class Hit implements Comparable { private final double adc; private final double time; - private Line3D wireLine; + private Line3D wireLine; private double phi; private double radius; private int nbOfWires; @@ -31,7 +30,7 @@ public class Hit implements Comparable { private double residual; private int trackId; - //updated constructor with ADC + //updated constructor with ADC public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca, double _ADC, double _Time) { this.id = _Id; this.superLayerId = _Super_layer; @@ -91,13 +90,13 @@ public double getDoca() { return doca; } - public Line3D getLine() { - return wireLine; - } - - public double getRadius() { - return radius; - } + public Line3D getLine() { + return wireLine; + } + + public double getRadius() { + return radius; + } public int getNbOfWires() { return nbOfWires; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java index 8bdddfeff7..a31775aa58 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/HitReader.java @@ -35,7 +35,7 @@ public final void fetch_AHDCHits(DataEvent event, AlertDCDetector detector) { } RawDataBank bankDGTZ = new RawDataBank("AHDC::adc"); - bankDGTZ.read(event); + bankDGTZ.read(event); for (int i = 0; i < bankDGTZ.rows(); i++) { int id = bankDGTZ.trueIndex(i) + 1; @@ -82,8 +82,8 @@ public final void fetch_AHDCHits(DataEvent event, AlertDCDetector detector) { double doca = p0 + p1*Math.pow(time,1.0) + p2*Math.pow(time,2.0) + p3*Math.pow(time,3.0) + p4*Math.pow(time,4.0) + p5*Math.pow(time, 5.0); if (time < 0) doca = 0; Hit h = new Hit(id, superlayer, layer, wire, doca, adc, time); - h.setWirePosition(detector); - hits.add(h); + h.setWirePosition(detector); + hits.add(h); } } } From 0ab72d1e619a7ac9ccb85133d93ce56afc87271a Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 13:44:48 +0100 Subject: [PATCH 15/34] implement a better comparable method for Hit --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 30 ++++++++++++++----- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 13 ++++++-- 2 files changed, 32 insertions(+), 11 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 e9e9d84115..647333a74a 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 @@ -1,13 +1,12 @@ package org.jlab.rec.ahdc.Hit; -import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; -import org.jlab.geom.prim.Line3D; -import org.jlab.geom.prim.Point3D; - 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.detector.alert.AHDC.AlertDCDetector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; public class Hit implements Comparable { @@ -61,12 +60,27 @@ public String toString() { return "Hit{" + "_Super_layer=" + superLayerId + ", _Layer=" + layerId + ", _Wire=" + wireId + ", _Doca=" + doca + ", _Phi=" + phi + '}'; } + // Should return + // 0 if equality + // +1 if this is bigger than arg0 + // -1 if this is lower than arg0 @Override public int compareTo(Hit arg0) { - if (this.superLayerId == arg0.superLayerId && this.layerId == arg0.layerId && this.wireId == arg0.wireId) { - return 0; - } else { - return 1; + if (this.superLayerId == arg0.superLayerId && this.layerId == arg0.layerId) { // same layer, so they have the same nbOfWires + if (this.wireId == 1 && arg0.wireId == this.nbOfWires) { + return 1; + } + else if (this.wireId == this.nbOfWires && arg0.wireId == 1) { + return -1; + } + else { + return Integer.compare(this.wireId, arg0.wireId); + } + } + else { + int this_value = 10*this.superLayerId + this.layerId; + int value = 10*arg0.superLayerId + arg0.layerId; + return Integer.compare(this_value, value); } } 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 e0dbb8f6dc..bea3eb96c5 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 @@ -14,6 +14,7 @@ import org.jlab.rec.ahdc.Track.Track; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.KFMonitor; +import java.util.Collections; import java.util.ArrayList; import java.util.HashMap; @@ -77,9 +78,15 @@ private void propagation(ArrayList tracks, DataEvent event, final double double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit ArrayList AHDC_hits = track.getHits(); - //----------------------------------------- - // add a routine to sort the list of hits - //----------------------------------------- + System.out.println("*** Before sorting ***"); + for (Hit hit : AHDC_hits) { + System.out.println(hit); + } + Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java + System.out.println("*** After sorting ***"); + for (Hit hit : AHDC_hits) { + System.out.println(hit); + } double zbeam = 0; if(IsVtxDefined)zbeam = vz_constraint;//test From 82d0d68c36258891d7c8233155be8fa39cad006d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:07:20 +0100 Subject: [PATCH 16/34] fix distance for the beamline --- .../main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 5 +++++ 1 file changed, 5 insertions(+) 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 4e5fc3be09..72c2a89760 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 @@ -3,6 +3,7 @@ import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; import org.jlab.rec.ahdc.Hit.Hit; // A weird object that just want to be considered as a Hit @@ -34,4 +35,8 @@ public RealVector get_Vector_beam() { public Line3D getLine() { return beamline; } + + public double distance(Point3D point3D) { + return this.beamline.distance(point3D).length(); + } } From 8d20ed500819086ce14b8cfae71fa10f805cce66 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:07:40 +0100 Subject: [PATCH 17/34] update maesurement noise --- .../alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 647333a74a..824271d3cc 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 @@ -169,7 +169,7 @@ public RealVector get_Vector() { } public RealMatrix get_MeasurementNoise() { - return new Array2DRowRealMatrix(new double[][]{{0.0225}}); + return new Array2DRowRealMatrix(new double[][]{{0.09}}); } public RealVector get_Vector_beam() { From a5093c842c0b24129e1004978fdbcba356fe5c20 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:08:24 +0100 Subject: [PATCH 18/34] fix residual and update Niter --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 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 bea3eb96c5..a27cf03589 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 @@ -35,7 +35,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 10; // number of iterations of the Kalman Filter + private final int Niter = 60; // number of iterations of the Kalman Filter private final boolean IsVtxDefined = false; // never used so far // can be useful to take into account the vertex of the electron @@ -78,15 +78,20 @@ private void propagation(ArrayList tracks, DataEvent event, final double double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Initialization hit ArrayList AHDC_hits = track.getHits(); - System.out.println("*** Before sorting ***"); + /*System.out.println("*** Before sorting ***"); for (Hit hit : AHDC_hits) { System.out.println(hit); - } + }*/ + /*System.out.println("*** Set Random position ***"); + Collections.shuffle(AHDC_hits); + for (Hit hit : AHDC_hits) { + System.out.println(hit); + }*/ Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java - System.out.println("*** After sorting ***"); + /*System.out.println("*** After sorting ***"); for (Hit hit : AHDC_hits) { System.out.println(hit); - } + }*/ double zbeam = 0; if(IsVtxDefined)zbeam = vz_constraint;//test @@ -113,7 +118,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); // Loop over number of iterations for (int k = 0; k < Niter; k++) { - //System.out.println("--------- ForWard propagation !! ---------"); + //System.out.println("Niter:" + k + " --------- ForWard propagation !! ---------"); //Reset error covariance: //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed for (Indicator indicator : forwardIndicators) { @@ -133,7 +138,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double } } - //System.out.println("--------- BackWard propagation !! ---------"); + //System.out.println("Niter:" + k + "--------- BackWard propagation !! ---------"); for (Indicator indicator : backwardIndicators) { TrackFitter.predict(indicator); // KFMonitor: save state and error covariance matrix @@ -158,7 +163,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if (indicator.haveAHit()) { if( indicator.hit.getId()>0){ - indicator.hit.setResidual(TrackFitter.residual(indicator)); + indicator.hit.setResidual(PostFitPropagator.residual(indicator)); } } } From e04f4ef8ee3939512a4bdd1dee6080405b10b258 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:08:57 +0100 Subject: [PATCH 19/34] tmp: remove banks before addign new ones --- .../alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java index 3dccfefdee..9a7c6c30b6 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java @@ -274,7 +274,8 @@ public int compare(Hit a1, Hit a2) { DataBank recoKFTracksBank = writer.fillAHDCKFTrackBank(event, AHDC_Tracks); DataBank AIPredictionBanks = writer.fillAIPrediction(event, predictions); DataBank recoKFMonitorBank = writer.fillKFMonitorBank(event, AHDC_Tracks); - + + event.removeBanks("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction","AHDC::kftrack:mon"); event.appendBank(recoHitsBank); event.appendBank(recoPreClusterBank); event.appendBank(recoClusterBank); From 43bf801161c54e7aa33f3d036f96a0978d5eabb4 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:50:00 +0100 Subject: [PATCH 20/34] fix hit beam id --- .../src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 72c2a89760..34b01638cb 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 @@ -19,7 +19,7 @@ public class Hit_beam extends Hit { // FOr clarity, we should remove them //public Hit_beam(double x, double y , double z) {} public Hit_beam(double x, double y , double z) { - super(-1, -1, -1, -1, 0, 0, -1); // just a line parallel to the beam axis, defined by two end points + super(0,0,0,0, 0, 0, -1); // just a line parallel to the beam axis, defined by two end points this.x = x; this.y = y; this.z = z; From b368e267b8372bf53b296945015aa0b0f3b47b40 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 15:50:24 +0100 Subject: [PATCH 21/34] fix indicator unique --- .../jlab/rec/ahdc/KalmanFilter/Indicator.java | 31 +++++++++---------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java index d4f88e55d7..1193676060 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java @@ -47,23 +47,20 @@ public boolean haveAHit() { // It is just for testing // Only used by KFMonitor public int getUniqueId() { - if (hit != null) { - return (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(); - } - else if (Math.abs(R) < 1e-9) { - // beamline - return 0; - } - else if (Math.abs(R - 3.0) < 1e-9) { - // inner face of the target straw - return 1; - } - else if (Math.abs(R - 3.060) < 1e-9) { - // outer face of the target straw - return 2; - } - else { - return -1; + if (hit != null) { + return (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(); + // for the beamline it will return 0 } + else if (Math.abs(R - 3.0) < 1e-9) { + // inner face of the target straw + return 1; + } + else if (Math.abs(R - 3.060) < 1e-9) { + // outer face of the target straw + return 2; + } + else { + return -1; + } } } From 491a7dba78e9825bff9889e52902d2638a66d6d2 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 16:48:19 +0100 Subject: [PATCH 22/34] add debug test for hit comparison --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) 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 824271d3cc..02c2e1b9b9 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 @@ -4,9 +4,11 @@ import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; +import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; public class Hit implements Comparable { @@ -180,4 +182,20 @@ public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); } + public static void main(String[] args) { + AlertDCDetector factory = (new AlertDCFactory()).createDetectorCLAS(new DatabaseConstantProvider()); + System.out.println("Run test: comparison between two hits."); + Hit h1 = new Hit(1,1,1,1,0,0,0); + Hit h2 = new Hit(1,1,1,47,0,0,0); + Hit h3 = new Hit(1,2,1,47,0,0,0); + h1.setWirePosition(factory); + h2.setWirePosition(factory); + System.out.println("h1 : " + h1); + System.out.println("h2 : " + h2); + System.out.println("numWires : " + h1.getNbOfWires()); + System.out.println("h1 compare to h2 : " + h1.compareTo(h2)); + System.out.println("h2 compare to h1 : " + h2.compareTo(h1)); + System.out.println("h1 compare to h3 : " + h1.compareTo(h3)); + } + } From fc49c2dbe7cb243f6ebc6b873c79be09dd86237c Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 24 Nov 2025 16:55:02 +0100 Subject: [PATCH 23/34] update test --- .../alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java | 1 + 1 file changed, 1 insertion(+) 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 02c2e1b9b9..8999d6580e 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 @@ -192,6 +192,7 @@ public static void main(String[] args) { h2.setWirePosition(factory); System.out.println("h1 : " + h1); System.out.println("h2 : " + h2); + System.out.println("h3 : " + h3); System.out.println("numWires : " + h1.getNbOfWires()); System.out.println("h1 compare to h2 : " + h1.compareTo(h2)); System.out.println("h2 compare to h1 : " + h2.compareTo(h1)); From ec63ad76d9aec1cad3a6a8347afc9e23fa0c6e82 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 25 Nov 2025 11:55:43 +0100 Subject: [PATCH 24/34] feat!: The Kalman Filter uses all available hits From 987e5fd926083d832f02a8393df403db62bc89e7 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 26 Nov 2025 13:24:10 +0100 Subject: [PATCH 25/34] use electron vertex --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 17 ++++++++++++++--- 1 file changed, 14 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 a27cf03589..1877340923 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 @@ -36,8 +36,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} private final int Niter = 60; // number of iterations of the Kalman Filter - private final boolean IsVtxDefined = false; // never used so far - // can be useful to take into account the vertex of the electron + private boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -51,6 +50,18 @@ private void propagation(ArrayList tracks, DataEvent event, final double final double[] B = {0.0, 0.0, magfield / 10 * tesla}; //DataBank mcBank = event.getBank("MC::Particle"); + // Load electron vertex + if (event.hasBank("REC::Particle")) { + DataBank recBank = event.getBank("REC::Particle"); + int row = 0; + while ((!IsVtxDefined) && row < recBank.rows()) { + if (recBank.getInt("pid", row) == 11) { + IsVtxDefined = true; + vz_constraint = recBank.getFloat("vz",row); + } + row++; + } + } // Initialization material map HashMap materialHashMap = materialGeneration(); @@ -94,7 +105,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double }*/ double zbeam = 0; - if(IsVtxDefined)zbeam = vz_constraint;//test + if(IsVtxDefined)zbeam = vz_constraint; // Define forward and backward indicator // cf. Indicator.java final ArrayList forwardIndicators = forwardIndicators(AHDC_hits, materialHashMap); From 7af331f68507ca488e8091f13d17dba00b4db46e Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 9 Dec 2025 13:56:16 +0100 Subject: [PATCH 26/34] do not use the electron vertex yet --- .../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 1877340923..9b58983c8e 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 @@ -51,7 +51,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double //DataBank mcBank = event.getBank("MC::Particle"); // Load electron vertex - if (event.hasBank("REC::Particle")) { + /*if (event.hasBank("REC::Particle")) { DataBank recBank = event.getBank("REC::Particle"); int row = 0; while ((!IsVtxDefined) && row < recBank.rows()) { @@ -61,7 +61,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double } row++; } - } + }*/ // Initialization material map HashMap materialHashMap = materialGeneration(); From 72b955c76660560f9b5d4bd2847ae263989db70b Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Thu, 11 Dec 2025 16:36:55 +0100 Subject: [PATCH 27/34] Fix the Kalmam Filter propagator - have a better control of the number of steps in the propagator - get rid of the notion of *Indicator* - propagate between hits - read materials from the MaterialMap class --- .../jlab/rec/ahdc/KalmanFilter/Indicator.java | 66 ----- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 112 +++---- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 215 ++++---------- .../rec/ahdc/KalmanFilter/Propagator.java | 279 ++++++++++-------- .../rec/ahdc/KalmanFilter/RungeKutta4.java | 1 - .../jlab/rec/ahdc/KalmanFilter/Stepper.java | 8 +- .../org/jlab/rec/ahdc/Track/KFMonitor.java | 2 +- 7 files changed, 266 insertions(+), 417 deletions(-) delete mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java deleted file mode 100644 index 1193676060..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Indicator.java +++ /dev/null @@ -1,66 +0,0 @@ -package org.jlab.rec.ahdc.KalmanFilter; - -import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.rec.ahdc.Hit.Hit; - -// Give the relative positions of the main constituents of the detector -// - beamline -// - target (inner face) -// - target (outer face) -// - ordered list of hits (wires) -// - hit(s) in layer 1 -// - hit(s) in layer 2 -// - -// - hit(s) in layer 8 -// It is used to indicate the Stepper(.java) -// - where to go -// - the medium/material it will cross -// - the step length to be used by RungeKutta4 (RK4) -public class Indicator { - - public double R; // radius position of the detector constituents (cylindrical geometry) - public double h; // step length used by RK4 - public Hit hit; // is null when the stepper moves toward target surfaces - public boolean direction; // true is the forward propagation (beamline to last layer) and false is the backward propagation (last layer to beamline) - public Material material; // medium property - - public Indicator(double r, double h, Hit hit, boolean direction, Material material) { - this.R = r; - this.h = h; - this.hit = hit; - this.direction = direction; - this.material = material; - } - - // Not well named - // Tell if we have a measurement or not at this indication (indicator) - // - Each AHDC hit (sense wire) provides a measurement, it is the distance (doca, measured distance to the wire) - // - When the stepper move toward the beamline (backward propagation), we also have a measurement (the distance to the beamline is 0) - public boolean haveAHit() { - boolean res = this.hit != null; - if (this.R == 0 && !direction) res = true; // beamline during backward propagation - return res; - } - - - // This method is not self consistent. We should not have hardcoded values - // It is just for testing - // Only used by KFMonitor - public int getUniqueId() { - if (hit != null) { - return (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(); - // for the beamline it will return 0 - } - else if (Math.abs(R - 3.0) < 1e-9) { - // inner face of the target straw - return 1; - } - else if (Math.abs(R - 3.060) < 1e-9) { - // outer face of the target straw - return 2; - } - else { - return -1; - } - } -} 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 aa682d4fea..01dd67457b 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 @@ -1,11 +1,15 @@ package org.jlab.rec.ahdc.KalmanFilter; +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; import org.jlab.clas.pdg.PhysicsConstants; +import org.jlab.clas.tracking.kalmanfilter.Material; import org.jlab.geom.prim.Point3D; +import org.jlab.rec.ahdc.Hit.Hit; public class KFitter { @@ -13,29 +17,32 @@ public class KFitter { private RealMatrix errorCovariance; public final Stepper stepper; private final Propagator propagator; + private final HashMap materialHashMap; public double chi2 = 0; // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; private boolean isvertexdefined = false; - public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator) { + public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator, final HashMap materialHashMap) { this.stateEstimation = initialStateEstimate; this.errorCovariance = initialErrorCovariance; this.stepper = stepper; this.propagator = propagator; + this.materialHashMap = materialHashMap; } - public void predict(Indicator indicator) throws Exception { + public void predict(Hit hit, boolean direction) throws Exception { // Initialization - stepper.initialize(indicator); + stepper.initialize(direction); Stepper stepper1 = new Stepper(stepper.y); + stepper1.initialize(direction); // project the state estimation ahead (a priori state) : xHat(k)- = f(xHat(k-1)) - stateEstimation = propagator.f(stepper, indicator); + stateEstimation = propagator.f(stepper, hit, materialHashMap); // project the covariance matrix ahead - RealMatrix transitionMatrix = F(indicator, stepper1); + RealMatrix transitionMatrix = F(hit, stepper1); RealMatrix transitionMatrixT = transitionMatrix.transpose(); double px = Math.abs(stepper.y[3]); @@ -58,7 +65,7 @@ public void predict(Indicator indicator) throws Exception { double dE = Math.abs(stepper.dEdx); double K = 0.000307075; - double sigma2_dE = indicator.material.getDensity() * K * indicator.material.getZoverA() / beta2 * tmax * s / 10 * (1.0 - beta2 / 2) * 1000 * 1000;//in MeV^2 + double sigma2_dE = stepper.material.getDensity() * K * stepper.material.getZoverA() / beta2 * tmax * s / 10 * (1.0 - beta2 / 2) * 1000 * 1000;//in MeV^2 double dp_prim_ddE = (E + dE) / Math.sqrt((E + dE) * (E + dE) - mass * mass); double sigma2_px = Math.pow(px / p, 2) * Math.pow(dp_prim_ddE, 2) * sigma2_dE; double sigma2_py = Math.pow(py / p, 2) * Math.pow(dp_prim_ddE, 2) * sigma2_dE; @@ -72,12 +79,13 @@ public void predict(Indicator indicator) throws Exception { errorCovariance = (transitionMatrix.multiply(errorCovariance.multiply(transitionMatrixT))).add(processNoise); } - public void correct(Indicator indicator) { + public void correct(Hit hit) { RealVector z; RealMatrix measurementNoise; RealMatrix measurementMatrix; RealVector h; - if (indicator.R == 0.0 && !indicator.direction) { + // check if the hit is the beamline + if (hit.getRadius() < 1) { double z_beam_res_sq = 1.e10;//in mm if(isvertexdefined)z_beam_res_sq = 4.0;//assuming 2. mm resolution measurementNoise = @@ -89,12 +97,14 @@ public void correct(Indicator indicator) { });//3x3 measurementMatrix = H_beam(stateEstimation);//6x3 h = h_beam(stateEstimation);//3x1 - z = indicator.hit.get_Vector_beam();//0! - } else { - measurementNoise = indicator.hit.get_MeasurementNoise();//1x1 - measurementMatrix = H(stateEstimation, indicator);//6x1 - h = h(stateEstimation, indicator);//1x1 - z = indicator.hit.get_Vector();//1x1 + z = hit.get_Vector_beam();//0! + } + // else, it is an AHDC hits + else { + measurementNoise = hit.get_MeasurementNoise();//1x1 + measurementMatrix = H(stateEstimation, hit);//6x1 + h = h(stateEstimation, hit);//1x1 + z = hit.get_Vector();//1x1 } RealMatrix measurementMatrixT = measurementMatrix.transpose(); @@ -122,41 +132,41 @@ public void correct(Indicator indicator) { stepper.y = stateEstimation.toArray(); } - public double residual(Indicator indicator) { - double d = indicator.hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); - return indicator.hit.getDoca()-d; + public double residual(Hit hit) { + double d = hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); + return hit.getDoca()-d; } public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } - private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { + private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { - double[] dfdx = subfunctionF(indicator, stepper1, 0); - double[] dfdy = subfunctionF(indicator, stepper1, 1); - double[] dfdz = subfunctionF(indicator, stepper1, 2); - double[] dfdpx = subfunctionF(indicator, stepper1, 3); - double[] dfdpy = subfunctionF(indicator, stepper1, 4); - double[] dfdpz = subfunctionF(indicator, stepper1, 5); + double[] dfdx = subfunctionF(hit, stepper1, 0); + double[] dfdy = subfunctionF(hit, stepper1, 1); + double[] dfdz = subfunctionF(hit, stepper1, 2); + double[] dfdpx = subfunctionF(hit, stepper1, 3); + double[] dfdpy = subfunctionF(hit, stepper1, 4); + double[] dfdpz = subfunctionF(hit, stepper1, 5); return MatrixUtils.createRealMatrix(new double[][]{ {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(Indicator indicator, Stepper stepper1, int i) throws Exception { + double[] subfunctionF(Hit 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); - stepper_plus.initialize(indicator); - stepper_minus.initialize(indicator); + stepper_plus.initialize(stepper1.direction); + stepper_minus.initialize(stepper1.direction); stepper_plus.y[i] = stepper_plus.y[i] + h; stepper_minus.y[i] = stepper_minus.y[i] - h; - propagator.f(stepper_plus, indicator); - propagator.f(stepper_minus, indicator); + propagator.f(stepper_plus, hit, materialHashMap); + propagator.f(stepper_minus, hit, materialHashMap); double dxdi = (stepper_plus.y[0] - stepper_minus.y[0]) / (2 * h); double dydi = (stepper_plus.y[1] - stepper_minus.y[1]) / (2 * h); @@ -168,28 +178,27 @@ private RealMatrix F(Indicator indicator, Stepper stepper1) throws Exception { return new double[]{dxdi, dydi, dzdi, dpxdi, dpydi, dpzdi}; } - //measurement matrix in 1 dimension: minimize distance - doca - private RealVector h(RealVector x, Indicator indicator) { - double d = indicator.hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + // Measurement matrix in 1x1 dimension: minimize distance - doca + private RealVector h(RealVector x, Hit hit) { + double d = hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); return MatrixUtils.createRealVector(new double[]{d}); } - //measurement matrix in 1 dimension: minimize distance - doca - private RealMatrix H(RealVector x, Indicator indicator) { + // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) + private RealMatrix H(RealVector x, Hit hit) { - double ddocadx = subfunctionH(x, indicator, 0); - double ddocady = subfunctionH(x, indicator, 1); - double ddocadz = subfunctionH(x, indicator, 2); - double ddocadpx = subfunctionH(x, indicator, 3); - double ddocadpy = subfunctionH(x, indicator, 4); - double ddocadpz = subfunctionH(x, indicator, 5); + 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); - // As per my understanding: ddocadx,y,z -> = dr/dx,y,z, etc return MatrixUtils.createRealMatrix(new double[][]{ {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); } - double subfunctionH(RealVector x, Indicator indicator, int i) { + double subfunctionH(RealVector x, Hit hit, int i) { double h = 1e-8;// in mm RealVector x_plus = x.copy(); RealVector x_minus = x.copy(); @@ -197,12 +206,13 @@ private RealMatrix H(RealVector x, Indicator indicator) { x_plus.setEntry(i, x_plus.getEntry(i) + h); x_minus.setEntry(i, x_minus.getEntry(i) - h); - double doca_plus = h(x_plus, indicator).getEntry(0); - double doca_minus = h(x_minus, indicator).getEntry(0); + 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); @@ -214,7 +224,8 @@ private RealVector h_beam(RealVector x) { 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); @@ -249,21 +260,12 @@ private RealMatrix H_beam(RealVector x) { }); } - /** - * Returns a copy of the current state estimation vector. - * - * @return the state estimation vector - */ public RealVector getStateEstimationVector() { return stateEstimation.copy(); } public RealMatrix getErrorCovarianceMatrix() { return errorCovariance.copy(); - } - - public double getMomentum() { - return Math.sqrt(stateEstimation.getEntry(3) * stateEstimation.getEntry(3) + stateEstimation.getEntry(4) * stateEstimation.getEntry(4) + stateEstimation.getEntry(5) * stateEstimation.getEntry(5)); } public void setVertexDefined(boolean isvtxdef) {isvertexdefined = isvtxdef;} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index 9b58983c8e..709a02e346 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 @@ -7,7 +7,6 @@ import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.clas.tracking.kalmanfilter.Units; import org.jlab.geom.prim.Point3D; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; @@ -19,23 +18,13 @@ import java.util.ArrayList; import java.util.HashMap; -/** - * TODO : - Fix multi hit on the same layer - * - Optimize measurement noise and probably use doca as weight - * - Fix the wire number (-1) - * - use a flag for simulation - * - add target in the map - * - move map to initialization engine - * - flag for target material - * - error px0 use MC !! Bad !! FIX IT FAST - */ // masses/energies should be in MeV; distances should be in mm public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 60; // number of iterations of the Kalman Filter + private final int Niter = 1; // number of iterations of the Kalman Filter private boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -48,6 +37,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; + HashMap materialHashMap = MaterialMap.generateMaterials(); //DataBank mcBank = event.getBank("MC::Particle"); // Load electron vertex @@ -62,16 +52,13 @@ private void propagation(ArrayList tracks, DataEvent event, final double row++; } }*/ - - // Initialization material map - HashMap materialHashMap = materialGeneration(); + // Loop over tracks - // Each track candidate form the HelixFitter is independently processed by the Kalman Filter int trackId = 0; for (Track track : tracks) { trackId++; track.set_trackId(trackId); - // Initialization State Vector + // Initialize state vector double x0 = 0.0; double y0 = 0.0; double z0 = track.get_Z0(); @@ -87,96 +74,80 @@ private void propagation(ArrayList tracks, DataEvent event, final double pz0 = mcBank.getFloat("pz", trackId-1)*1000; }*/ double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; - // Initialization hit + // Read list of hits ArrayList AHDC_hits = track.getHits(); - /*System.out.println("*** Before sorting ***"); - for (Hit hit : AHDC_hits) { - System.out.println(hit); - }*/ - /*System.out.println("*** Set Random position ***"); - Collections.shuffle(AHDC_hits); - for (Hit hit : AHDC_hits) { - System.out.println(hit); - }*/ Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java - /*System.out.println("*** After sorting ***"); - for (Hit hit : AHDC_hits) { - System.out.println(hit); - }*/ double zbeam = 0; if(IsVtxDefined)zbeam = vz_constraint; - // Define forward and backward indicator - // cf. Indicator.java - final ArrayList forwardIndicators = forwardIndicators(AHDC_hits, materialHashMap); - final ArrayList backwardIndicators = backwardIndicators(AHDC_hits, materialHashMap, zbeam); // Start propagation Stepper stepper = new Stepper(y); RungeKutta4 RK4 = new RungeKutta4(proton, 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); - //first 3 lines in mm^2; last 3 lines in MeV^2 RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); - KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator); + KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); TrackFitter.setVertexDefined(IsVtxDefined); - // KFmonitor: save initial state and error covariance matrix + // KFmonitor: save state and error covariance matrix track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - // Loop over number of iterations + + // Loop over number of iterations for (int k = 0; k < Niter; k++) { - //System.out.println("Niter:" + k + " --------- ForWard propagation !! ---------"); //Reset error covariance: //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed - for (Indicator indicator : forwardIndicators) { - // Prediction - TrackFitter.predict(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - // I don't see the utility of this - if( k==0 && indicator.hit.getId()>0){ // first iteration and indicator != beamline (because its id is -1) - indicator.hit.setResidualPrefit(TrackFitter.residual(indicator)); - } - // Correction only if we have a measure (hit) - TrackFitter.correct(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, k, 0, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - } + + // Forward propagation + //System.out.println("==================================> Forward Progation"); + for (Hit hit : AHDC_hits) { + TrackFitter.predict(hit, true); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + // I don't see the utility of this + // so I comment it, we will also need to delete the residual_prefit attribut in the AHDC::hits bank + /*if( k==0 && indicator.hit.getId()>0){ // first iteration and indicator != beamline (because its id is -1) + indicator.hit.setResidualPrefit(TrackFitter.residual(indicator)); + }*/ + TrackFitter.correct(hit); + track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + } + //System.out.println("==================================> Backward Progation"); + // Backward propagation (last layer to first layer) + for (int i = AHDC_hits.size() - 2; i >= 0; i--) { + Hit hit = AHDC_hits.get(i); + TrackFitter.predict(hit, false); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + TrackFitter.correct(hit); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + } + // Backward propagation (first layer to beamline) + { + Hit hit = new Hit_beam(0, 0, zbeam); + TrackFitter.predict(hit, false); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + TrackFitter.correct(hit); + track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + } - //System.out.println("Niter:" + k + "--------- BackWard propagation !! ---------"); - for (Indicator indicator : backwardIndicators) { - TrackFitter.predict(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - TrackFitter.correct(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, k, 1, indicator.getUniqueId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - } - } + } RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumForKF(x_out); - + //System.out.println("==================================> PostFit Progation"); //Residual, path and AHDC exit momentum calculation post fit: - KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4)); - for (Indicator indicator : forwardIndicators) { - PostFitPropagator.predict(indicator); - // KFMonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, indicator.getUniqueId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - if (indicator.haveAHit()) { - if( indicator.hit.getId()>0){ - indicator.hit.setResidual(PostFitPropagator.residual(indicator)); - } - } + KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); + for (Hit hit : AHDC_hits) { + PostFitPropagator.predict(hit, true); + track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); + if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits + hit.setResidual(PostFitPropagator.residual(hit)); + } } // Fill track and hit bank @@ -201,92 +172,8 @@ private void propagation(ArrayList tracks, DataEvent event, final double }//end of loop on track candidates } catch (Exception e) { // e.printStackTrace(); + System.out.println("===> Error"); } } - private HashMap materialGeneration() { - Units units = Units.CM; - - String name_De = "deuteriumGas"; - double thickness_De = 1; - double density_De = 9.37E-4;// 5.5 atm - double ZoverA_De = 0.496499; - double X0_De = 1.3445E+5; // I guess X0 is not even used??? - double IeV_De = 19.2; - - org.jlab.clas.tracking.kalmanfilter.Material deuteriumGas = new org.jlab.clas.tracking.kalmanfilter.Material(name_De, thickness_De, density_De, ZoverA_De, X0_De, IeV_De, units); - - String name_Bo = "BONuS12Gas";//80% He, 20% CO2 - double thickness_Bo = 1; - double density_Bo = 1.39735E-3; - double ZoverA_Bo = 0.49983; - double X0_Bo = 3.69401E+4; - double IeV_Bo = 73.5338; - - org.jlab.clas.tracking.kalmanfilter.Material BONuS12 = new org.jlab.clas.tracking.kalmanfilter.Material(name_Bo, thickness_Bo, density_Bo, ZoverA_Bo, X0_Bo, IeV_Bo, units); - - String name_My = "Mylar"; - double thickness_My = 1; - double density_My = 1.4; - double ZoverA_My = 0.52037; - double X0_My = 28.54; - double IeV_My = 78.7; - - org.jlab.clas.tracking.kalmanfilter.Material Mylar = new org.jlab.clas.tracking.kalmanfilter.Material(name_My, thickness_My, density_My, ZoverA_My, X0_My, IeV_My, units); - - String name_Ka = "Kapton"; - double thickness_Ka = 1; - double density_Ka = 1.42; - double ZoverA_Ka = 0.51264; - double X0_Ka = 28.57; - double IeV_Ka = 79.6; - - org.jlab.clas.tracking.kalmanfilter.Material Kapton = new org.jlab.clas.tracking.kalmanfilter.Material(name_Ka, thickness_Ka, density_Ka, ZoverA_Ka, X0_Ka, IeV_Ka, units); - - return new HashMap() { - { - put("deuteriumGas", deuteriumGas); - put("Kapton", Kapton); - put("Mylar", Mylar); - put("BONuS12Gas", BONuS12); - } - }; - } - - ArrayList forwardIndicators(ArrayList hitArrayList, HashMap materialHashMap) { - ArrayList forwardIndicators = new ArrayList<>(); - //R, h, defined in mm! - forwardIndicators.add(new Indicator(3.0, 0.2, null, true, materialHashMap.get("deuteriumGas"))); - forwardIndicators.add(new Indicator(3.060, 0.001, null, true, materialHashMap.get("Kapton"))); - for (Hit hit : hitArrayList) { - forwardIndicators.add(new Indicator(hit.getRadius(), 0.1, hit, true, materialHashMap.get("BONuS12Gas"))); - } - return forwardIndicators; - } - - ArrayList backwardIndicators(ArrayList hitArrayList, HashMap materialHashMap) { - ArrayList backwardIndicators = new ArrayList<>(); - //R, h, defined in mm! - for (int i = hitArrayList.size() - 2; i >= 0; i--) { - backwardIndicators.add(new Indicator(hitArrayList.get(i).getRadius(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); - } - backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); - backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); - Hit hit = new Hit_beam(0, 0, 0); - backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); - return backwardIndicators; - } - - ArrayList backwardIndicators(ArrayList hitArrayList, HashMap materialHashMap, double vz) { - ArrayList backwardIndicators = new ArrayList<>(); - //R, h, defined in mm! - for (int i = hitArrayList.size() - 2; i >= 0; i--) { - backwardIndicators.add(new Indicator(hitArrayList.get(i).getRadius(), 0.1, hitArrayList.get(i), false, materialHashMap.get("BONuS12Gas"))); - } - backwardIndicators.add(new Indicator(3.060, 1, null, false, materialHashMap.get("BONuS12Gas"))); - backwardIndicators.add(new Indicator(3.0, 0.001, null, false, materialHashMap.get("Kapton"))); - Hit hit = new Hit_beam(0, 0, vz); - backwardIndicators.add(new Indicator(0.0, 0.2, hit, false, materialHashMap.get("deuteriumGas"))); - return backwardIndicators; - } } 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 fe03799437..ef70d4b7b9 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,16 +1,15 @@ 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 java.io.FileWriter; -import java.io.IOException; -import java.io.Writer; -import java.util.Arrays; +import org.jlab.clas.tracking.kalmanfilter.Material; +import java.util.HashMap; +import org.jlab.rec.ahdc.Hit.Hit; // All distances here should be in mm. -// Do all those hardcoded values even make sense??? public class Propagator { private final RungeKutta4 RK4; @@ -19,145 +18,175 @@ public Propagator(RungeKutta4 rungeKutta4) { this.RK4 = rungeKutta4; } - // Propagate the stepper toward the next indicator - void propagate(Stepper stepper, Indicator indicator) { - // ------------------------------------------------------------ - final int maxNbOfStep = 10000; // do not allow more than 10000 steps (very critical cases) - final double R = indicator.R; + // This function needs to know: + // - the direction + // These information are already accessible in the stepper, but where shoulb specifiy them + // Propagate the stepper toward the next hit + //void propagate(Stepper stepper, Indicator indicator, HashMap materialHashMap) { + void propagate(Stepper stepper, Hit hit, HashMap materialHashMap) { + + // Do not allow more than 10000 steps (very critical cases) + final int maxNbOfStep = 10000; + + // Initialize a stepper, used to save the previous stepper + Stepper prevStepper = new Stepper(stepper.y); + // Initialize distances, d should always disminish double dMin = Double.MAX_VALUE; - double d = 0; + double prev_dMin = Double.MAX_VALUE; + double d; + + // Intialize radius + double R = stepper.r(); + double prev_R = R; + + // Initialize the stepper size + stepper.h = 0.1; + + // Initialize material + // R = 3 mm is the location of the target + if (R < 3) { + stepper.material = materialHashMap.get("deuteriumGas"); + } else { + stepper.material = materialHashMap.get("BONuS12Gas"); + } - // System.out.println("R = " + R); - // stepper.print(); - - // We should use a "while" instead of "for" - for (int nbStep = 0; nbStep < maxNbOfStep; nbStep++) { - double previous_r = stepper.r(); + // boolean : check if we reach or not the target + boolean target_reached = false; + boolean target_crossed = false; + double thickness = 0; + + // Do the propagation + int nbStep = 0; + while (nbStep < maxNbOfStep) { + nbStep++; + // Save previous state + prevStepper.y = Arrays.copyOf(stepper.y, stepper.y.length); + prev_R = stepper.r(); + // Take a step RK4.doOneStep(stepper); - double r = stepper.r(); - // stepper.print(); - - - if (stepper.direction) { // forward propagation - if (r >= R - 2 - 1.5 * stepper.h) stepper.h = 1e-2; // ? - if (indicator.hit != null) { // the indicator is a hit/wire - if (r >= R - 2) { // only when the stepper is 2 mm closer to the layer of the hit - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); // distance of the stepper with respect to the wire - if (d < dMin) dMin = d; // the distance d should continue to diminish - } - if (r >= R + 2 || d > dMin) { // r should not be bigger than R - // if the distance d starts to grow again, stop the propagation! - // in that case, the good stepper is the previous one! - // System.out.println("dMin = " + dMin); - break; - } - } else { // the indicator is the target face (inner or outer) - if (r >= R) { - break; // break as soon as r >= R, again, may the good stepper is the previous one + // compute distance with respect to the hit + // this distance should always disminish + d = hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); + R = stepper.r(); + // check the evolution + if (d < dMin) { + prev_dMin = dMin; + dMin = d; + } + else { + // go back to the previous step and stop the propagation + stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + // --------------------------------------- + // this part can be optomized + // we can make the step size dynamical + // e.g : if d >= dMin, go back to the previous step and reduce the step size by a factor 5 + // to be done later + // --------------------------------------- + break; + } + + // When the propagation is done between the beamline and the first AHDC hit, + // we should take care of the target material + // target is located at R = 3 mm with a thickness of 0.060 mm (60 µm) + // if R < 5 mm, be careful! + if (R < 5) { + + // We try to cross the target + if (((prev_R < 3 && R > 3) || (prev_R > 3 && R < 3)) && !target_reached) { + // We want to cross the target with a very small step < 0.060 mm + if (Math.abs(R - prev_R) > 0.060) { + stepper.h /= 5; + // redo the propagation + // so, go back to the previous step + stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + dMin = prev_dMin; + continue; } - } - } else { // backward direction - if (r <= R + 2 + 1.5 * stepper.h) stepper.h = 1e-2; // ? - if (R < 5) { // we are now close the target (outer face is at 3.060 mm) - if (stepper.h < 1e-4) stepper.h = 1e-4; // this step is too short // it increases the computing time - if ((previous_r - r) < 0) { // in the backward propagation, we expect r to decrease; if not stop the propagation! - break; + // we can consider that we have reached the target and that we are ready to cross it + else { + target_reached = true; + // redo the propagation + // so, go back to the previous step + stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + dMin = prev_dMin; + // We are now ready to cross the target + // set the target material + stepper.material = materialHashMap.get("Kapton"); + // update the step size + stepper.h = 0.003; + // initialise the crossed thickness + thickness = 0; + continue; } } - if (indicator.hit != null) { // wire/hit or beamline - if (r <= R + 2) { // here R < r and r decreases, we enter here when r becomes 2 mm closer to R - // we compute the distance, we expect it to decrease - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); - if (d < dMin) dMin = d; - } - if (r <= R - 2 || d > dMin) { // if the distance grows again, stop the propagation - // System.out.println("dMin = " + dMin); - break; - } - } else { // the indicator is the target face (inner or outer) - if (r <= R) { - break; // break as soon as r <= R, again, may the good stepper is the previous one + + // cross the target for 0.06 mm + if (target_reached && !target_crossed) { + thickness += Math.abs(R - prev_R); + if (thickness > 0.06) { + target_crossed = true; + // update stepper size + // we go back to a normal propagation + stepper.h = 0.1; + // update stepper material + if (stepper.direction) { // forward propagation + stepper.material = materialHashMap.get("BONuS12Gas"); + } + else { // backward propagation + stepper.material = materialHashMap.get("deuteriumGas"); + } } } - - + } - - } + //System.out.printf("nbstep : %d (%s)\n", nbStep, stepper.direction ? "forward" : "backward"); + } - public RealVector f(Stepper stepper, Indicator indicator) { - propagate(stepper, indicator); + public RealVector f(Stepper stepper, Hit hit, HashMap materialHashMap) { + propagate(stepper, hit, materialHashMap); return new ArrayRealVector(stepper.y); } - - // ------------------------------- - // not used, to be deleted - // ------------------------------- - public void propagateAndWrite(Stepper stepper, Indicator indicator, Writer writer) { - // ------------------------------------------------------------ - final int maxNbOfStep = 10000; + +} + +// Draft +// ------------------------------------------------------------ + /*final int maxNbOfStep = 10000; // do not allow more than 10000 steps (very critical cases) final double R = indicator.R; double dMin = Double.MAX_VALUE; double d = 0; + Stepper prevStepper = new Stepper(stepper.y); + prevStepper.initialize(indicator); - System.out.println("R = " + R); - stepper.print(); - try {writer.write("" + Arrays.toString(stepper.y) + '\n');} catch (Exception e) {e.printStackTrace();} - - for (int nbStep = 0; nbStep < maxNbOfStep; nbStep++) { - double previous_r = stepper.r(); + int nbStep = 0; + while (nbStep < maxNbOfStep) { + nbStep++; + prevStepper.y = Arrays.copyOf(stepper.y, stepper.y.length); RK4.doOneStep(stepper); - double r = stepper.r(); - stepper.print(); - try {writer.write("" + Arrays.toString(stepper.y) + '\n');} catch (Exception e) {e.printStackTrace();} - - if (stepper.direction) { - if (r >= R - 2 - 1.5 * stepper.h) stepper.h = 1e-2; - if (indicator.hit != null) { - if (r >= R - 2) { - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); - System.out.println("d = " + d); - if (d < dMin) dMin = d; - } - if (r >= R + 2 || d > dMin) { - System.out.println("dMin = " + dMin); - break; - } - } else { - if (r >= R) { - break; - } - } - } else { - if (r <= R + 2 + 1.5 * stepper.h) stepper.h = 1e-2; - if (R < 5) { - if (stepper.h < 1e-4) stepper.h = 1e-4; - if ((previous_r - r) < 0) { - break; - } + + // Compute the distance + if (indicator.hit != null) { // the indicator is a hit/wire + d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); // distance of the stepper with respect to the wire + } else { // the indicator is the target face (inner or outer) + double d0 = Math.abs(indicator.R - stepper.r()); + d = Math.max(d0, stepper.h); + if (flag) { + System.out.printf("[%4d] r : %f ---> R : %2.5f (%s)\n", nbStep, stepper.r(), indicator.R, stepper.direction ? "forward" : "backward"); } - if (indicator.hit != null) { - if (r <= R + 2) { - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); - if (d < dMin) dMin = d; - } - if (r <= R - 2 || d > dMin) { - System.out.println("dMin = " + dMin); - break; - } - } else { - if (r <= R) { - break; - } - } - - - } - } - } -} + } + + // The first time, it will disminish because dMin = "Infinity" + // Starting nbStep 2, the distance should continue to disminish. If not, stop the propagation + if (d < dMin) { + dMin = d; + } else { + stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + break; + } + + }*/ \ No newline at end of file diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java index e66f1a3a56..410217ebc4 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java @@ -3,7 +3,6 @@ import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.pdg.PhysicsConstants; -import java.util.Arrays; // Compute the next trajectory point using the equation of motion // the path length is given by the stepper diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java index 1d6b52709c..fa78ade7fd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java @@ -23,14 +23,12 @@ public Stepper(double[] y) { // The initialisation tells the stepper where to go (i.e toward the `indicator`) // Cf. notes in Indicator.java // The initilisation is done before any propagation between indicators, in the predict() method of KFitter.java - public void initialize(Indicator indicator) { + //public void initialize(Indicator indicator) { + public void initialize(boolean direction) { this.s = 0; this.dEdx = 0; - this.h = indicator.h; - this.material = indicator.material; - this.direction = indicator.direction; + this.direction = direction; this.is_in_drift = false; - } // radius, distance to the beamline diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java index f96cb33d4a..93c5d488fd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java @@ -7,7 +7,7 @@ public class KFMonitor { private int trackid; //< trackid private int Niter; //< iteration number of the Kalman Filter algorithm private int orientation; //< forward (0), backward (1), postfit (2) propagation - private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 ; target straw surface ==> 1 (inner face), 2 (outer face) [[ remark: if indicator < 10 then it is not a wire ]] + private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 private int status; //< state just after: prediction (0) or correction (1) private RealVector state; //< vector containing x, y, z, px, py, pz private RealMatrix errorCovarianceMatrix; // error covariance matrix, only the diagonal matter (x-x')^2, (y-y')^2, ..., (pz-pz')^2 From 94e578e2249df6bb7f41a7c5bfba60799e2eb682 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 15 Dec 2025 11:05:42 +0100 Subject: [PATCH 28/34] uncomment error status --- .../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 709a02e346..7a51ecb055 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 @@ -171,8 +171,8 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.set_n_hits(AHDC_hits.size()); }//end of loop on track candidates } catch (Exception e) { - // e.printStackTrace(); - System.out.println("===> Error"); + e.printStackTrace(); + System.out.println("======> Kalman Filter Error"); } } From 7f7a6255e926af613a633ddc9f084cf9af09008c Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Mon, 15 Dec 2025 11:05:58 +0100 Subject: [PATCH 29/34] try to optimize the propagator --- .../rec/ahdc/KalmanFilter/Propagator.java | 113 ++++++++---------- .../jlab/rec/ahdc/KalmanFilter/Stepper.java | 34 ++++++ 2 files changed, 83 insertions(+), 64 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 ef70d4b7b9..7e188ca92b 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,6 +1,6 @@ package org.jlab.rec.ahdc.KalmanFilter; -import java.util.Arrays; +//import java.util.Arrays; import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; @@ -18,22 +18,20 @@ public Propagator(RungeKutta4 rungeKutta4) { this.RK4 = rungeKutta4; } - // This function needs to know: - // - the direction - // These information are already accessible in the stepper, but where shoulb specifiy them // Propagate the stepper toward the next hit - //void propagate(Stepper stepper, Indicator indicator, HashMap materialHashMap) { void propagate(Stepper stepper, Hit hit, HashMap materialHashMap) { // Do not allow more than 10000 steps (very critical cases) final int maxNbOfStep = 10000; - // Initialize a stepper, used to save the previous stepper - Stepper prevStepper = new Stepper(stepper.y); + // Initialise steppers to save the two previous states + Stepper prev_stepper = new Stepper(stepper); + Stepper prev_prev_stepper = new Stepper(stepper); // Initialize distances, d should always disminish double dMin = Double.MAX_VALUE; double prev_dMin = Double.MAX_VALUE; + double prev_prev_dMin = Double.MAX_VALUE; double d; // Intialize radius @@ -41,7 +39,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM double prev_R = R; // Initialize the stepper size - stepper.h = 0.1; + stepper.h = 0.5; // Initialize material // R = 3 mm is the location of the target @@ -57,33 +55,55 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM double thickness = 0; // Do the propagation + //stepper.print(); int nbStep = 0; while (nbStep < maxNbOfStep) { nbStep++; // Save previous state - prevStepper.y = Arrays.copyOf(stepper.y, stepper.y.length); + prev_prev_stepper.copyContent(prev_stepper); + prev_stepper.copyContent(stepper); prev_R = stepper.r(); // Take a step RK4.doOneStep(stepper); // compute distance with respect to the hit // this distance should always disminish d = hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); + /*if (hit.getRadius() < 34 && hit.getRadius() > 30 && !stepper.direction) { // beamline for backward propagation + System.out.print("d = " + d + " "); + stepper.print(); + }*/ R = stepper.r(); // check the evolution if (d < dMin) { + prev_prev_dMin = prev_dMin; prev_dMin = dMin; dMin = d; } else { - // go back to the previous step and stop the propagation - stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); - // --------------------------------------- - // this part can be optomized - // we can make the step size dynamical - // e.g : if d >= dMin, go back to the previous step and reduce the step size by a factor 5 - // to be done later - // --------------------------------------- - break; + // reduce the step size + //if (stepper.h > 0.05) { + //if (Math.abs(prev_dMin - dMin) > 0.05) { + if (stepper.distance(prev_stepper) > 0.05) { + // we didn't reach the precision, so + // go back two steps before + // and redo the propagation + stepper.copyContent(prev_prev_stepper); + dMin = prev_prev_dMin; + stepper.h /= 2; + continue; + } + else { + // we reach the precision, so + // go back to the previous step and stop the propagation + stepper.copyContent(prev_stepper); + stepper.h = 0.5; + break; + } + // we reach the precision, so + // go back to the previous step and break + /*stepper.copyContent(prev_stepper); + stepper.h = 0.1; + break;*/ } // When the propagation is done between the beamline and the first AHDC hit, @@ -96,10 +116,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM if (((prev_R < 3 && R > 3) || (prev_R > 3 && R < 3)) && !target_reached) { // We want to cross the target with a very small step < 0.060 mm if (Math.abs(R - prev_R) > 0.060) { - stepper.h /= 5; // redo the propagation // so, go back to the previous step - stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + stepper.copyContent(prev_stepper); + stepper.h /= 5; dMin = prev_dMin; continue; } @@ -108,12 +128,12 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM target_reached = true; // redo the propagation // so, go back to the previous step - stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); + stepper.copyContent(prev_stepper); dMin = prev_dMin; // We are now ready to cross the target // set the target material stepper.material = materialHashMap.get("Kapton"); - // update the step size + // copyContent the step size stepper.h = 0.003; // initialise the crossed thickness thickness = 0; @@ -126,10 +146,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM thickness += Math.abs(R - prev_R); if (thickness > 0.06) { target_crossed = true; - // update stepper size + // copyContent stepper size // we go back to a normal propagation - stepper.h = 0.1; - // update stepper material + stepper.h = 0.5; + // copyContent stepper material if (stepper.direction) { // forward propagation stepper.material = materialHashMap.get("BONuS12Gas"); } @@ -142,7 +162,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } } - //System.out.printf("nbstep : %d (%s)\n", nbStep, stepper.direction ? "forward" : "backward"); + /*System.out.printf("nbstep : %d (%s)\n", nbStep, stepper.direction ? "forward" : "backward"); + if (hit.getId() == 0) { + System.out.printf("******** Beamline <-- radius = %f, stepper.h = %f\n", stepper.r(), stepper.h); + }*/ } @@ -152,41 +175,3 @@ public RealVector f(Stepper stepper, Hit hit, HashMap material } } - -// Draft -// ------------------------------------------------------------ - /*final int maxNbOfStep = 10000; // do not allow more than 10000 steps (very critical cases) - final double R = indicator.R; - - double dMin = Double.MAX_VALUE; - double d = 0; - Stepper prevStepper = new Stepper(stepper.y); - prevStepper.initialize(indicator); - - int nbStep = 0; - while (nbStep < maxNbOfStep) { - nbStep++; - prevStepper.y = Arrays.copyOf(stepper.y, stepper.y.length); - RK4.doOneStep(stepper); - - // Compute the distance - if (indicator.hit != null) { // the indicator is a hit/wire - d = indicator.hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); // distance of the stepper with respect to the wire - } else { // the indicator is the target face (inner or outer) - double d0 = Math.abs(indicator.R - stepper.r()); - d = Math.max(d0, stepper.h); - if (flag) { - System.out.printf("[%4d] r : %f ---> R : %2.5f (%s)\n", nbStep, stepper.r(), indicator.R, stepper.direction ? "forward" : "backward"); - } - } - - // The first time, it will disminish because dMin = "Infinity" - // Starting nbStep 2, the distance should continue to disminish. If not, stop the propagation - if (d < dMin) { - dMin = d; - } else { - stepper.y = Arrays.copyOf(prevStepper.y, prevStepper.y.length); - break; - } - - }*/ \ No newline at end of file diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java index fa78ade7fd..3b99e41182 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java @@ -19,6 +19,32 @@ public class Stepper { public Stepper(double[] y) { this.y = Arrays.copyOf(y, y.length); } + + // Copy constructor + public Stepper(Stepper stepper) { + this.y = Arrays.copyOf(stepper.y, stepper.y.length); + this.h = stepper.h; + this.material = stepper.material; + this.s = stepper.s; + this.sTot = stepper.sTot; + this.s_drift = stepper.s_drift; + this.is_in_drift = stepper.is_in_drift; + this.dEdx = stepper.dEdx; + this.direction = stepper.direction; + } + + // Copy the contents but do not change the reference + public void copyContent(Stepper stepper) { + this.y = Arrays.copyOf(stepper.y, stepper.y.length); + this.h = stepper.h; + this.material = stepper.material; + this.s = stepper.s; + this.sTot = stepper.sTot; + this.s_drift = stepper.s_drift; + this.is_in_drift = stepper.is_in_drift; + this.dEdx = stepper.dEdx; + this.direction = stepper.direction; + } // The initialisation tells the stepper where to go (i.e toward the `indicator`) // Cf. notes in Indicator.java @@ -49,4 +75,12 @@ public void print() { public String toString() { return "" + sTot + ' ' + y[0] + ' ' + y[1] + ' ' + y[2] + ' ' + y[3] + ' ' + y[4] + ' ' + y[5]; } + + public double distance(Stepper stepper) { + double dx = this.y[0] - stepper.y[0]; + double dy = this.y[1] - stepper.y[1]; + double dz = this.y[2] - stepper.y[2]; + return Math.sqrt(dx*dx + dy*dy + dz*dz); + } + } From f3cdf3193be857be6cfa0f91eb59b6900788f419 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 16 Dec 2025 11:52:49 +0100 Subject: [PATCH 30/34] increase initial covariance errors --- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 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 7a51ecb055..fc1646d08e 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,5 +1,9 @@ package org.jlab.rec.ahdc.KalmanFilter; +import java.util.ArrayList; +import java.util.Collections; +import java.util.HashMap; + import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; @@ -7,16 +11,12 @@ import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.geom.prim.Point3D; -import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; -import org.jlab.rec.ahdc.Track.Track; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Track.KFMonitor; -import java.util.Collections; +import org.jlab.rec.ahdc.Track.Track; -import java.util.ArrayList; -import java.util.HashMap; +//import org.apache.commons.math3.linear.RealMatrixFormat; // masses/energies should be in MeV; distances should be in mm @@ -24,7 +24,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 1; // number of iterations of the Kalman Filter + private final int Niter = 60; // number of iterations of the Kalman Filter private boolean IsVtxDefined = false; private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -89,8 +89,8 @@ private void propagation(ArrayList tracks, DataEvent event, final double // 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[][]{{1.00, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 1.00, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 25.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 25.0}}); - KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); + RealMatrix initialErrorCovariance = MatrixUtils.createRealMatrix(new double[][]{{50.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 50.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 900.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 100.00, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 100.00, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 900.0}}); + KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); TrackFitter.setVertexDefined(IsVtxDefined); // KFmonitor: save state and error covariance matrix @@ -98,24 +98,14 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Loop over number of iterations for (int k = 0; k < Niter; k++) { - //Reset error covariance: - //TrackFitter.ResetErrorCovariance(initialErrorCovariance); // that can be very interesting, to be checked later, this has an effect on the convergence speed - // Forward propagation - //System.out.println("==================================> Forward Progation"); for (Hit hit : AHDC_hits) { TrackFitter.predict(hit, true); track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - // I don't see the utility of this - // so I comment it, we will also need to delete the residual_prefit attribut in the AHDC::hits bank - /*if( k==0 && indicator.hit.getId()>0){ // first iteration and indicator != beamline (because its id is -1) - indicator.hit.setResidualPrefit(TrackFitter.residual(indicator)); - }*/ TrackFitter.correct(hit); track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } - //System.out.println("==================================> Backward Progation"); // Backward propagation (last layer to first layer) for (int i = AHDC_hits.size() - 2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); @@ -133,9 +123,19 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } - } + /*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.setPositionAndMomentumForKF(x_out); @@ -171,8 +171,8 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.set_n_hits(AHDC_hits.size()); }//end of loop on track candidates } catch (Exception e) { - e.printStackTrace(); - System.out.println("======> Kalman Filter Error"); + //e.printStackTrace(); + //System.out.println("======> Kalman Filter Error"); } } From 155a51306525cd65fffa026e1d300afb33751ffc Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 16 Dec 2025 12:26:10 +0100 Subject: [PATCH 31/34] remove KFMonitor --- etc/bankdefs/util/bankSplit.py | 2 +- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 42 ------------------- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 36 +--------------- .../org/jlab/rec/ahdc/Track/KFMonitor.java | 33 --------------- .../java/org/jlab/rec/ahdc/Track/Track.java | 7 ---- .../org/jlab/service/ahdc/AHDCEngine.java | 6 +-- 6 files changed, 5 insertions(+), 121 deletions(-) delete mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java diff --git a/etc/bankdefs/util/bankSplit.py b/etc/bankdefs/util/bankSplit.py index 74aedb3cab..9340d248f5 100755 --- a/etc/bankdefs/util/bankSplit.py +++ b/etc/bankdefs/util/bankSplit.py @@ -65,7 +65,7 @@ def create(dirname, banklist): raster = ["RASTER::position"] rich = ["RICH::tdc","RICH::Ring","RICH::Particle"] rtpc = ["RTPC::hits","RTPC::tracks","RTPC::KFtracks"] -alert = ["AHDC::track", "AHDC::mc", "AHDC::hits", "AHDC::preclusters", "AHDC::clusters", "AHDC::kftrack", "AHDC::ai:prediction", "AHDC::kftrack:mon"] +alert = ["AHDC::track", "AHDC::mc", "AHDC::hits", "AHDC::preclusters", "AHDC::clusters", "AHDC::kftrack", "AHDC::ai:prediction"] dets = band + raster + rich + rtpc + alert # additions for the calibration schema: 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 a7885abd3d..fcf11e9222 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,7 +7,6 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.Track.Track; -import org.jlab.rec.ahdc.Track.KFMonitor; import org.apache.commons.math3.linear.RealVector; import org.apache.commons.math3.linear.RealMatrix; @@ -183,45 +182,4 @@ public DataBank fillAIPrediction(DataEvent event, ArrayList pre return bank; } - public DataBank fillKFMonitorBank(DataEvent event, ArrayList tracks) { - - int nentries = 0; - for (Track track : tracks) { - nentries += track.get_ListOfKFMonitors().size(); - } - - DataBank bank = event.createBank("AHDC::kftrack:mon", nentries); - - int row = 0; - - for (Track track : tracks) { - if (track == null) continue; - for (KFMonitor monitor : track.get_ListOfKFMonitors()) { - RealVector state = monitor.get_state(); - RealMatrix errorCovarianceMatrix = monitor.get_errorCovarianceMatrix(); - bank.setFloat("x", row, (float) state.getEntry(0)); - bank.setFloat("y", row, (float) state.getEntry(1)); - bank.setFloat("z", row, (float) state.getEntry(2)); - bank.setFloat("px", row, (float) state.getEntry(3)); - bank.setFloat("py", row, (float) state.getEntry(4)); - bank.setFloat("pz", row, (float) state.getEntry(5)); - bank.setShort("trackid", row, (short) monitor.get_trackid()); - bank.setShort("Niter", row, (short) monitor.get_Niter()); - bank.setShort("orientation", row, (short) monitor.get_orientation()); - bank.setShort("indicator", row, (short) monitor.get_indicator()); - bank.setShort("status", row, (short) monitor.get_status()); - bank.setFloat("var_x", row, (float) errorCovarianceMatrix.getEntry(0,0)); - bank.setFloat("var_y", row, (float) errorCovarianceMatrix.getEntry(1,1)); - bank.setFloat("var_z", row, (float) errorCovarianceMatrix.getEntry(2,2)); - bank.setFloat("var_px", row, (float) errorCovarianceMatrix.getEntry(3,3)); - bank.setFloat("var_py", row, (float) errorCovarianceMatrix.getEntry(4,4)); - bank.setFloat("var_pz", row, (float) errorCovarianceMatrix.getEntry(5,5)); - - row++; - } - } - - return bank; - } - } 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 fc1646d08e..bd478e7166 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 @@ -13,7 +13,6 @@ import org.jlab.clas.tracking.kalmanfilter.Material; import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; -import org.jlab.rec.ahdc.Track.KFMonitor; import org.jlab.rec.ahdc.Track.Track; //import org.apache.commons.math3.linear.RealMatrixFormat; @@ -38,21 +37,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; HashMap materialHashMap = MaterialMap.generateMaterials(); - - //DataBank mcBank = event.getBank("MC::Particle"); - // Load electron vertex - /*if (event.hasBank("REC::Particle")) { - DataBank recBank = event.getBank("REC::Particle"); - int row = 0; - while ((!IsVtxDefined) && row < recBank.rows()) { - if (recBank.getInt("pid", row) == 11) { - IsVtxDefined = true; - vz_constraint = recBank.getFloat("vz",row); - } - row++; - } - }*/ - + // Loop over tracks int trackId = 0; for (Track track : tracks) { @@ -65,14 +50,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double double px0 = track.get_px(); double py0 = track.get_py(); double pz0 = track.get_pz(); - // using or not mc - // will be deleted in the final code - /*if (IsMC) { - z0 = mcBank.getFloat("vz", trackId-1)*10; - px0 = mcBank.getFloat("px", trackId-1)*1000; - py0 = mcBank.getFloat("py", trackId-1)*1000; - pz0 = mcBank.getFloat("pz", trackId-1)*1000; - }*/ + double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; // Read list of hits ArrayList AHDC_hits = track.getHits(); @@ -93,34 +71,25 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, stepper, propagator, materialHashMap); TrackFitter.setVertexDefined(IsVtxDefined); - // KFmonitor: save state and error covariance matrix - track.add_KFMonitor(new KFMonitor(trackId, 0, 0, 0, 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); - // Loop over number of iterations for (int k = 0; k < Niter; k++) { // Forward propagation for (Hit hit : AHDC_hits) { TrackFitter.predict(hit, true); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); TrackFitter.correct(hit); - track.add_KFMonitor(new KFMonitor(trackId, k, 0, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } // Backward propagation (last layer to first layer) for (int i = AHDC_hits.size() - 2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); TrackFitter.predict(hit, false); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); TrackFitter.correct(hit); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } // Backward propagation (first layer to beamline) { Hit hit = new Hit_beam(0, 0, zbeam); TrackFitter.predict(hit, false); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); TrackFitter.correct(hit); - track.add_KFMonitor(new KFMonitor(trackId, k, 1, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 1, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); } } @@ -144,7 +113,6 @@ private void propagation(ArrayList tracks, DataEvent event, final double KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); - track.add_KFMonitor(new KFMonitor(trackId, Niter, 2, (hit.getSuperLayerId()*10 + hit.getLayerId())*100 + hit.getWireId(), 0, TrackFitter.getStateEstimationVector(), TrackFitter.getErrorCovarianceMatrix())); if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits hit.setResidual(PostFitPropagator.residual(hit)); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java deleted file mode 100644 index 93c5d488fd..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/KFMonitor.java +++ /dev/null @@ -1,33 +0,0 @@ -package org.jlab.rec.ahdc.Track; - -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.linear.RealMatrix; - -public class KFMonitor { - private int trackid; //< trackid - private int Niter; //< iteration number of the Kalman Filter algorithm - private int orientation; //< forward (0), backward (1), postfit (2) propagation - private int indicator; //< wire ==> layer*100 + component ; beamline ==> 0 - private int status; //< state just after: prediction (0) or correction (1) - private RealVector state; //< vector containing x, y, z, px, py, pz - private RealMatrix errorCovarianceMatrix; // error covariance matrix, only the diagonal matter (x-x')^2, (y-y')^2, ..., (pz-pz')^2 - - // constructor - public KFMonitor(int _trackid, int _Niter, int _orientation, int _indicator, int _status, RealVector _state, RealMatrix _errorCovarianceMatrix) { - this.trackid = _trackid; - this.Niter = _Niter; - this.orientation = _orientation; - this.indicator = _indicator; - this.status = _status; - this.state = _state.copy(); - this.errorCovarianceMatrix = _errorCovarianceMatrix.copy(); - } - - public int get_trackid() {return trackid;} - public int get_Niter() {return Niter;} - public int get_orientation() {return orientation;} - public int get_indicator() {return indicator;} - public int get_status() {return status;} - public RealVector get_state() {return state.copy();} - public RealMatrix get_errorCovarianceMatrix() {return errorCovarianceMatrix.copy();} -} 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 279ccf0fd9..ea1548a083 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 @@ -41,9 +41,6 @@ public class Track { private double dEdx_kf = 0; ///< deposited energy per path length (adc/mm) private double p_drift_kf = 0; ///< momentum in the drift region (MeV) private double path_kf = 0; ///< length of the track (mm) - // Monitoring of the Kalman Filter - private List ListOfKFMonitors = new ArrayList<>(); - public Track(List clusters) { this._Clusters = clusters; @@ -204,8 +201,4 @@ public double getPz0_kf() { public double get_p_drift_kf() {return p_drift_kf;} public double get_path_kf() {return path_kf;} - // KF monitoring - public void add_KFMonitor(KFMonitor _monitor) { ListOfKFMonitors.add(_monitor);} - public List get_ListOfKFMonitors() { return ListOfKFMonitors;} - } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java index 9a7c6c30b6..89c77c1ab0 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java @@ -111,7 +111,7 @@ public boolean init() { this.getConstantsManager().setVariation("default"); - this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction", "AHDC::kftrack:mon"); + this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction"); return true; } @@ -273,16 +273,14 @@ public int compare(Hit a1, Hit a2) { DataBank recoTracksBank = writer.fillAHDCTrackBank(event, AHDC_Tracks); DataBank recoKFTracksBank = writer.fillAHDCKFTrackBank(event, AHDC_Tracks); DataBank AIPredictionBanks = writer.fillAIPrediction(event, predictions); - DataBank recoKFMonitorBank = writer.fillKFMonitorBank(event, AHDC_Tracks); - event.removeBanks("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction","AHDC::kftrack:mon"); + event.removeBanks("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction"); event.appendBank(recoHitsBank); event.appendBank(recoPreClusterBank); event.appendBank(recoClusterBank); event.appendBank(recoTracksBank); event.appendBank(recoKFTracksBank); event.appendBank(AIPredictionBanks); - event.appendBank(recoKFMonitorBank); if (simulation) { DataBank recoMCBank = writer.fillAHDCMCTrackBank(event); From 48dc5b077363c8f59d1372f991944f220f050d30 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 16 Dec 2025 15:28:51 +0100 Subject: [PATCH 32/34] Update comments --- .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 1 + .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 28 +++++++---- .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 7 +++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 20 +++++--- .../rec/ahdc/KalmanFilter/MaterialMap.java | 13 ++--- .../rec/ahdc/KalmanFilter/Propagator.java | 48 ++++++++----------- .../rec/ahdc/KalmanFilter/RungeKutta4.java | 11 +++-- .../jlab/rec/ahdc/KalmanFilter/Stepper.java | 32 +++++++------ 8 files changed, 94 insertions(+), 66 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 8999d6580e..7e4f4546f5 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 @@ -174,6 +174,7 @@ public RealMatrix get_MeasurementNoise() { return new Array2DRowRealMatrix(new double[][]{{0.09}}); } + // a signature for KalmanFilter.Hit_beam public RealVector get_Vector_beam() { return null; } 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 34b01638cb..fd25c2e94d 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,37 +6,47 @@ import org.jlab.geom.prim.Point3D; import org.jlab.rec.ahdc.Hit.Hit; -// A weird object that just want to be considered as a Hit -// but it is not! -// It is due to the notion of Indicator that want to unify wire hits, beamline or target surfaces +/** + * A weird object that just want to be considered as a Hit + * The methods that matters are : distance(), getLine(), get_Vector_beam(), getRadius() + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class Hit_beam extends Hit { private double x,y,z; private double r,phi; Line3D beamline; - // A lot of argument are useless - // FOr clarity, we should remove them - //public Hit_beam(double x, double y , double z) {} public Hit_beam(double x, double y , double z) { - super(0,0,0,0, 0, 0, -1); // just a line parallel to the beam axis, defined by two end points + super(0,0,0,0, 0, 0, -1); this.x = x; this.y = y; this.z = z; this.r = Math.hypot(x,y); this.phi = Math.atan2(y,x); - beamline = new Line3D(x,y,0,x,y,1); + beamline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis } + @Override public RealVector get_Vector_beam() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } - + + @Override public Line3D getLine() { return beamline; } + @Override public double distance(Point3D point3D) { return this.beamline.distance(point3D).length(); } + + @Override + public double getRadius() { + return r; + } } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 01dd67457b..3f19dc0305 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 @@ -11,6 +11,13 @@ import org.jlab.geom.prim.Point3D; import org.jlab.rec.ahdc.Hit.Hit; +/** + * Implement the prediction and the correction stages of the Kalman Filter + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class KFitter { private RealVector stateEstimation; 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 bd478e7166..82866491cd 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 @@ -17,14 +17,21 @@ //import org.apache.commons.math3.linear.RealMatrixFormat; -// masses/energies should be in MeV; distances should be in mm - +/** + * This is the main routine of the Kalman Filter. The fit is done by a KFitter + * + * masses/energies should be in MeV; distances should be in mm + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 60; // number of iterations of the Kalman Filter - private boolean IsVtxDefined = false; + private final int Niter = 60; // number of iterations for the Kalman Filter + private boolean IsVtxDefined = false; // implemented but not used yet private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { @@ -108,8 +115,8 @@ private void propagation(ArrayList tracks, DataEvent event, final double RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumForKF(x_out); - //System.out.println("==================================> PostFit Progation"); - //Residual, path and AHDC exit momentum calculation post fit: + + // 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); for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); @@ -119,6 +126,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double } // 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; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java index fcfa011552..8dfd6eb1fb 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/MaterialMap.java @@ -5,12 +5,13 @@ import java.util.HashMap; -// In KalmanFilter.java, we have a refinition of this material -// We should merge these codes in the sense that KalmanFilter -// should refer to this class -// -// Also this is not exactly the material for ALERT -// To be modified +/** + * TO DO: use/rename the relevant materials for ALERT + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class MaterialMap { public static HashMap generateMaterials() { 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 7e188ca92b..e85021addb 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 @@ -9,7 +9,15 @@ import java.util.HashMap; import org.jlab.rec.ahdc.Hit.Hit; -// All distances here should be in mm. +/** + * Is responsible of the propagation + * + * All distances here should be in mm. + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class Propagator { private final RungeKutta4 RK4; @@ -55,11 +63,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM double thickness = 0; // Do the propagation - //stepper.print(); int nbStep = 0; while (nbStep < maxNbOfStep) { nbStep++; - // Save previous state + // Save previous states prev_prev_stepper.copyContent(prev_stepper); prev_stepper.copyContent(stepper); prev_R = stepper.r(); @@ -68,10 +75,6 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM // compute distance with respect to the hit // this distance should always disminish d = hit.distance(new Point3D(stepper.y[0], stepper.y[1], stepper.y[2])); - /*if (hit.getRadius() < 34 && hit.getRadius() > 30 && !stepper.direction) { // beamline for backward propagation - System.out.print("d = " + d + " "); - stepper.print(); - }*/ R = stepper.r(); // check the evolution if (d < dMin) { @@ -80,12 +83,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM dMin = d; } else { - // reduce the step size - //if (stepper.h > 0.05) { - //if (Math.abs(prev_dMin - dMin) > 0.05) { + // check the distance between the steps, they should be very big if (stepper.distance(prev_stepper) > 0.05) { - // we didn't reach the precision, so - // go back two steps before + // we didn't reach the precision + // so, go back two steps before, reduce the step size // and redo the propagation stepper.copyContent(prev_prev_stepper); dMin = prev_prev_dMin; @@ -93,22 +94,17 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM continue; } else { - // we reach the precision, so - // go back to the previous step and stop the propagation + // 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 stepper.copyContent(prev_stepper); stepper.h = 0.5; break; } - // we reach the precision, so - // go back to the previous step and break - /*stepper.copyContent(prev_stepper); - stepper.h = 0.1; - break;*/ } // When the propagation is done between the beamline and the first AHDC hit, // we should take care of the target material - // target is located at R = 3 mm with a thickness of 0.060 mm (60 µm) + // the target is located at R = 3 mm with a thickness of 0.060 mm (60 µm) // if R < 5 mm, be careful! if (R < 5) { @@ -117,7 +113,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM // We want to cross the target with a very small step < 0.060 mm if (Math.abs(R - prev_R) > 0.060) { // redo the propagation - // so, go back to the previous step + // so, go back to the previous step and reduce the step size stepper.copyContent(prev_stepper); stepper.h /= 5; dMin = prev_dMin; @@ -133,7 +129,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM // We are now ready to cross the target // set the target material stepper.material = materialHashMap.get("Kapton"); - // copyContent the step size + // update the step size stepper.h = 0.003; // initialise the crossed thickness thickness = 0; @@ -146,10 +142,10 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM thickness += Math.abs(R - prev_R); if (thickness > 0.06) { target_crossed = true; - // copyContent stepper size + // reset stepper size // we go back to a normal propagation stepper.h = 0.5; - // copyContent stepper material + // update stepper material if (stepper.direction) { // forward propagation stepper.material = materialHashMap.get("BONuS12Gas"); } @@ -162,10 +158,6 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } } - /*System.out.printf("nbstep : %d (%s)\n", nbStep, stepper.direction ? "forward" : "backward"); - if (hit.getId() == 0) { - System.out.printf("******** Beamline <-- radius = %f, stepper.h = %f\n", stepper.r(), stepper.h); - }*/ } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java index 410217ebc4..b192ca978e 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java @@ -4,9 +4,14 @@ import org.jlab.clas.pdg.PhysicsConstants; -// Compute the next trajectory point using the equation of motion -// the path length is given by the stepper -// doneOneStep() is the main method, it is used in the Propagator.propagate() method +/** + * Compute the next trajectory point using the equation of motion + * the step size is given by the stepper + * doneOneStep() is the main method, it is used in Propagator.propagate() + * + * @author Mathieu Ouillon + * @author Éric Fuchey + */ public class RungeKutta4 { private final int numberOfVariables; private final PDGParticle particle; diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java index 3b99e41182..1d0f6f63d8 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Stepper.java @@ -4,17 +4,23 @@ import java.util.Arrays; -// A trajectory point computed by the Kalman Filter +/** + * It is a screenshot of the system during the Kalman Filter procedure + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ public class Stepper { public double[] y; // position and momentum coordinates - public double h; // - public Material material; - public double s; // path length during the propagation between the current indicator and the next one given in initialize() - public double sTot = 0; // total lenght with respect to the beamline // increase in forward propagation, reach the maximum in the last layer and decrease in the backward propagation // cf doOneStep() method in RungeKutta4.java - public double s_drift = 0; // cf doOneStep() method in RungeKutta4.java - public boolean is_in_drift = false; // cf doOneStep() method in RungeKutta4.java - public double dEdx; // deposited energy during the propagation between the current indicator and the next one given in initialize() - public boolean direction; // true is the forward propagation (beamline to last layer) and false is the backward propagation (last layer to beamline) + public double h; // step size, managed in Propagator + public Material material; // managed in Propagator + public double s; + public double sTot = 0; + public double s_drift = 0; + public boolean is_in_drift = false; + public double dEdx; // s, sTot, s_drift, is_in_drift, dEdx are used in RungeKutta.doOneStep() + public boolean direction; // true (forward propagation), false (backward propagation) public Stepper(double[] y) { this.y = Arrays.copyOf(y, y.length); @@ -46,10 +52,7 @@ public void copyContent(Stepper stepper) { this.direction = stepper.direction; } - // The initialisation tells the stepper where to go (i.e toward the `indicator`) - // Cf. notes in Indicator.java - // The initilisation is done before any propagation between indicators, in the predict() method of KFitter.java - //public void initialize(Indicator indicator) { + // the initialisation is done in the KFitter.predict() public void initialize(boolean direction) { this.s = 0; this.dEdx = 0; @@ -57,7 +60,7 @@ public void initialize(boolean direction) { this.is_in_drift = false; } - // radius, distance to the beamline + // radius public double r() { return Math.hypot(y[0], y[1]); } @@ -76,6 +79,7 @@ public String toString() { return "" + sTot + ' ' + y[0] + ' ' + y[1] + ' ' + y[2] + ' ' + y[3] + ' ' + y[4] + ' ' + y[5]; } + // distance between tow steppers public double distance(Stepper stepper) { double dx = this.y[0] - stepper.y[0]; double dy = this.y[1] - stepper.y[1]; From 045081a6d3211aecc027efddd8a435a705a3047d Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 16 Dec 2025 09:46:14 -0500 Subject: [PATCH 33/34] reduce number of iterations --- .../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 82866491cd..805a736558 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 @@ -30,7 +30,7 @@ public class KalmanFilter { public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 60; // number of iterations for the Kalman Filter + private final int Niter = 40; // number of iterations for the Kalman Filter private boolean IsVtxDefined = false; // implemented but not used yet private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { From 1545c85ec7ee1728bd9e47da58d584df0f8f8948 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 16 Dec 2025 17:22:46 +0100 Subject: [PATCH 34/34] Reset submodule etc/data/magfield to development version --- etc/data/magfield | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/etc/data/magfield b/etc/data/magfield index a0e899ea96..ba858a8d51 160000 --- a/etc/data/magfield +++ b/etc/data/magfield @@ -1 +1 @@ -Subproject commit a0e899ea96c86609ea36f5111ed7c319e8463673 +Subproject commit ba858a8d51cc5bd0b49ce881b32d44aa57e3e68b