From 2f31211280435e3654677c79c54721747d21baa0 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 9 Jun 2026 09:27:51 -0400 Subject: [PATCH 1/8] extends the bank AHDC::kftrack and allow re-cooking of AHDC only --- etc/bankdefs/hipo4/alert.json | 68 +++++++++++++++++++ .../rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- .../org/jlab/service/ahdc/AHDCEngine.java | 1 + .../org/jlab/service/alert/ALERTEngine.java | 1 + 4 files changed, 71 insertions(+), 1 deletion(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 64d4faeec1..580a54aecc 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -408,6 +408,74 @@ "name": "sum_residuals", "type": "F", "info": "Sum of residuals (mm)" + }, { + "name": "atof_region", + "type": "I", + "info": "a flag to know if a track reach S1, S2, or S3; is n if Sn is reached; is 0 otherwise" + }, { + "name": "atof_s1_x", + "type": "F", + "info": "x position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_y", + "type": "F", + "info": "y position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_z", + "type": "F", + "info": "z position of the track if it reaches the surface s1" + }, { + "name": "atof_s1_path", + "type": "F", + "info": "total path from the beamline to the surface s1" + }, { + "name": "atof_s1_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s1" + }, { + "name": "atof_s2_x", + "type": "F", + "info": "x position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_y", + "type": "F", + "info": "y position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_z", + "type": "F", + "info": "z position of the track if it reaches the surface s2" + }, { + "name": "atof_s2_path", + "type": "F", + "info": "total path from the beamline to the surface s2" + }, { + "name": "atof_s2_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s2" + }, { + "name": "atof_s3_x", + "type": "F", + "info": "x position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_y", + "type": "F", + "info": "y position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_z", + "type": "F", + "info": "z position of the track if it reaches the surface s3" + }, { + "name": "atof_s3_path", + "type": "F", + "info": "total path from the beamline to the surface s3" + }, { + "name": "atof_s3_p", + "type": "F", + "info": "momentum whith which the track reaches the surface s3" + }, { + "name": "atof_match", + "type": "B", + "info": "1 if there is a ATOF match, 0 if not" } ] }, { 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 26dcf2ee7f..268e43862b 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 @@ -38,7 +38,7 @@ public class KalmanFilter { private boolean IsVtxDefined = false; // implemented but not used yet private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration - private double clas_alignement = -54; + private double clas_alignement = -75; public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { 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 b52d306f7a..715646a48a 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 @@ -98,6 +98,7 @@ public boolean processDataEventUser(DataEvent event) { ArrayList AHDC_Hits = hitReader.get_AHDCHits(); RecoBankWriter writer = new RecoBankWriter(); + event.removeBank("AHDC::hits"); DataBank recoHitsBank = writer.fillAHDCHitsBank(event, AHDC_Hits); if (recoHitsBank != null) event.appendBank(recoHitsBank); } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index c54a7636e5..9b2ba0f773 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -582,6 +582,7 @@ public boolean processDataEventUser(DataEvent event) { ///////////////////////////////////////////// // write the AHDC::kftrack bank in the event ///////////////////////////////////////////// + event.removeBank("AHDC::kftrack"); org.jlab.rec.ahdc.Banks.RecoBankWriter ahdc_writer = new org.jlab.rec.ahdc.Banks.RecoBankWriter(); DataBank recoKFTracksBank = ahdc_writer.fillAHDCKFTrackBank(event, AHDC_tracks); event.appendBank(recoKFTracksBank); From 257af6c56ef9117ca39352443ab454ac61f94f2f Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 9 Jun 2026 11:43:26 -0400 Subject: [PATCH 2/8] convenient use of PID in RK4 --- .../java/org/jlab/rec/ahdc/KalmanFilter/RungeKutta4.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 b192ca978e..0331dc0f93 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 @@ -103,8 +103,7 @@ public void doOneStep(Stepper stepper) { } private double[] f(double[] y) { - double charge = 1.0; // the charge should by given by the particle object - // for now, we assume it is a proton + double charge = particle.charge(); 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; @@ -133,13 +132,14 @@ private double[] f(double[] y) { private void energyLoss( double[] yIn, double h, org.jlab.clas.tracking.kalmanfilter.Material material) { double mass = particle.mass() * 1000; //particle mass defined in GeV, converted to MeV + int charge = particle.charge(); h /= 10; // h defined in mm, converted to cm double mom = Math.sqrt(yIn[3] * yIn[3] + yIn[4] * yIn[4] + yIn[5] * yIn[5]); double E = Math.sqrt(mom * mom + mass * mass); //material::getEloss(double p, double m) uses GeV and cm //see common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Material.java - double dedx = material.getEloss(mom/1000, mass/1000) * 1000;//Momentum, mass input in GeV, output in GeV/cm, converted to MeV/cm + double dedx = material.getEloss(mom/1000, mass/1000, charge) * 1000;//Momentum, mass input in GeV, output in GeV/cm, converted to MeV/cm double DeltaE = dedx * h; stepper.dEdx += DeltaE; From d8175bbdb3531c0bd0a731b095a6568f95d2c3ac Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 9 Jun 2026 13:11:19 -0400 Subject: [PATCH 3/8] create KFHit and refactor the kalman filter --- .../jlab/rec/ahdc/Banks/RecoBankWriter.java | 34 ++ .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 64 ++- .../jlab/rec/ahdc/KalmanFilter/Hit_beam.java | 52 --- .../org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 22 + .../jlab/rec/ahdc/KalmanFilter/KFitter.java | 140 +----- .../rec/ahdc/KalmanFilter/KalmanFilter.java | 411 ++++++++++++++---- .../rec/ahdc/KalmanFilter/Propagator.java | 10 +- .../rec/ahdc/KalmanFilter/RadialKFHit.java | 105 +++++ .../ahdc/KalmanFilter/RadialSurfaceKFHit.java | 79 ++++ .../java/org/jlab/rec/alert/Track/Track.java | 67 +++ .../org/jlab/service/alert/ALERTEngine.java | 196 ++++++++- 11 files changed, 890 insertions(+), 290 deletions(-) delete mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java create mode 100644 reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java index dead9dacdd..b276e4b605 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java @@ -7,6 +7,7 @@ import org.jlab.rec.ahdc.AHDCCluster.AHDCCluster; import org.jlab.rec.ahdc.DocaCluster.DocaCluster; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.alert.Track.Track; @@ -153,6 +154,39 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { bank.setFloat("chi2", row, (float) track.get_chi2()); bank.setFloat("sum_residuals", row, (float) track.get_sum_residuals()); + // track projection on ATOF surface S1 + Stepper stepper_s1 = track.get_ATOF_S1_stepper(); + if (stepper_s1 != null) { + bank.setFloat("atof_s1_x", row, (float) stepper_s1.y[0]); + bank.setFloat("atof_s1_y", row, (float) stepper_s1.y[1]); + bank.setFloat("atof_s1_z", row, (float) stepper_s1.y[2]); + bank.setFloat("atof_s1_path", row, (float) stepper_s1.sTot); + bank.setFloat("atof_s1_p", row, (float) stepper_s1.p()); + } + + // track projection on ATOF surface S2 + Stepper stepper_s2 = track.get_ATOF_S2_stepper(); + if (stepper_s2 != null) { + bank.setFloat("atof_s2_x", row, (float) stepper_s2.y[0]); + bank.setFloat("atof_s2_y", row, (float) stepper_s2.y[1]); + bank.setFloat("atof_s2_z", row, (float) stepper_s2.y[2]); + bank.setFloat("atof_s2_path", row, (float) stepper_s2.sTot); + bank.setFloat("atof_s2_p", row, (float) stepper_s2.p()); + } + + // track projection on ATOF surface S3 + Stepper stepper_s3 = track.get_ATOF_S3_stepper(); + if (stepper_s3 != null) { + bank.setFloat("atof_s3_x", row, (float) stepper_s3.y[0]); + bank.setFloat("atof_s3_y", row, (float) stepper_s3.y[1]); + bank.setFloat("atof_s3_z", row, (float) stepper_s3.y[2]); + bank.setFloat("atof_s3_path", row, (float) stepper_s3.sTot); + bank.setFloat("atof_s3_p", row, (float) stepper_s3.p()); + } + + bank.setInt("atof_region", row, track.get_ATOF_region()); + bank.setByte("atof_match", row, (byte) ((track.getATOFHits().size() > 0) ? 1 : 0)); + row++; } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java index 008e892104..18dd1bcad5 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Hit/Hit.java @@ -2,6 +2,7 @@ import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; import org.jlab.detector.calib.utils.DatabaseConstantProvider; @@ -9,8 +10,9 @@ import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; +import org.jlab.rec.ahdc.KalmanFilter.KFHit; -public class Hit implements Comparable { +public class Hit implements Comparable, KFHit { private final int id; private final int superLayerId; @@ -33,6 +35,16 @@ public class Hit implements Comparable { private int trackId; //updated constructor with ADC + /** + * + * @param _Id is AHDC::adc row id + 1 + * @param _Super_layer super layer id + * @param _Layer layer id + * @param _Wire wire id + * @param _Doca distance from the timing information using the time2distance, matches {@link Hit.#time} + * @param _ADC raw ADC + * @param _Time calibrated time + */ public Hit(int _Id, int _Super_layer, int _Layer, int _Wire, double _Doca, double _ADC, double _Time) { this.id = _Id; this.superLayerId = _Super_layer; @@ -110,6 +122,7 @@ public Line3D getLine() { return wireLine; } + @Override public double getRadius() { return radius; } @@ -136,6 +149,7 @@ public double getY() { public double getPhi() {return phi;} + /** Get calibrated ADC */ public double getADC() {return adc;} public double getResidual() { @@ -146,18 +160,22 @@ public void setResidual(double resid) { this.residual = resid; } + /** Set calibrated ToT */ public void setToT(double _tot) { this.tot = _tot; } + /** Get calibrated ToT */ public double getToT() { return tot; } + /** Set calibrated ADC */ public void setADC(double _adc) { this.adc = _adc; } + /** Get raw ADC */ public double getRawADC() { return raw_adc; } @@ -174,11 +192,13 @@ public void setTrackId(int _trackId) { this.trackId = _trackId; } - public RealVector get_Vector() { + @Override + public RealVector MeasurementVector() { return new ArrayRealVector(new double[]{this.doca}); } - public RealMatrix get_MeasurementNoise() { + @Override + public RealMatrix MeasurementNoiseMatrix() { double mean_error = 0.471; // mm (no difference between adc and time) double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm double error_on_time = (0.4423*time + 13.7215)/(0.846038*time + 31.9867); // mm @@ -188,11 +208,43 @@ public RealMatrix get_MeasurementNoise() { //return new Array2DRowRealMatrix(new double[][]{{0.09}}); } - // a signature for KalmanFilter.Hit_beam - public RealVector get_Vector_beam() { - return null; + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + double d = this.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); + return MatrixUtils.createRealVector(new double[]{d}); + } + + // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double ddocadx = partialProjectionMatrix(x, 0); + double ddocady = partialProjectionMatrix(x, 1); + double ddocadz = partialProjectionMatrix(x, 2); + double ddocadpx = partialProjectionMatrix(x, 3); + double ddocadpy = partialProjectionMatrix(x, 4); + double ddocadpz = partialProjectionMatrix(x, 5); + + return MatrixUtils.createRealMatrix(new double[][]{ + {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); + } + + private double partialProjectionMatrix(RealVector x, int i) { + double h = 1e-8;// in mm + RealVector x_plus = x.copy(); + RealVector x_minus = x.copy(); + + x_plus.setEntry(i, x_plus.getEntry(i) + h); + x_minus.setEntry(i, x_minus.getEntry(i) - h); + + double doca_plus = this.ProjectionFunction(x_plus).getEntry(0); + double doca_minus = this.ProjectionFunction(x_minus).getEntry(0); + + return (doca_plus - doca_minus) / (2 * h); } + @Override public double distance(Point3D point3D) { return this.wireLine.distance(point3D).length(); } 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 deleted file mode 100644 index fd25c2e94d..0000000000 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Hit_beam.java +++ /dev/null @@ -1,52 +0,0 @@ -package org.jlab.rec.ahdc.KalmanFilter; - -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.RealVector; -import org.jlab.geom.prim.Line3D; -import org.jlab.geom.prim.Point3D; -import org.jlab.rec.ahdc.Hit.Hit; - -/** - * 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; - - public Hit_beam(double x, double y , double z) { - super(0,0,0,0, 0, 0, -1); - this.x = x; - this.y = y; - this.z = z; - this.r = Math.hypot(x,y); - this.phi = Math.atan2(y,x); - beamline = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis - } - - @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/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java new file mode 100644 index 0000000000..3a2e693192 --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -0,0 +1,22 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Point3D; +/** + * An interface to unify the hits used the Kalman Filter (e.g AHDC hits, ATOF hits, beamline) + * + * @author Felix Touchte Codjo + */ +public interface KFHit { + public double distance(Point3D point3D); + public double getRadius(); + /** Return the measurement encoded in this KFHit */ + public RealVector MeasurementVector(); + /** Return the measurement noise matrix for this this KFHit */ + public RealMatrix MeasurementNoiseMatrix(); + /** Compute the measure for a given state vector */ + public RealVector ProjectionFunction(RealVector x); + /** Compute the Jacobian matrix of the {@link #ProjectionFunction(RealVector)} with respect of the components of the state vector */ + public RealMatrix ProjectionMatrix(RealVector x); +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 7263d6638f..8fb0ecaafd 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java @@ -2,7 +2,6 @@ import java.util.HashMap; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; @@ -22,24 +21,23 @@ public class KFitter { private RealVector stateEstimation; private RealMatrix errorCovariance; - public final Stepper stepper; + private Stepper stepper; private final Propagator propagator; private final HashMap materialHashMap; public double chi2 = 0; // masses/energies in MeV private final double electron_mass_c2 = PhysicsConstants.massElectron() * 1000; private final double proton_mass_c2 = PhysicsConstants.massProton() * 1000; - private double[] vertex_resolutions = {0.09,1e10}; // default values // dr^2 and dz^2 in mm^2 - public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Stepper stepper, final Propagator propagator, final HashMap materialHashMap) { - this.stateEstimation = initialStateEstimate; - this.errorCovariance = initialErrorCovariance; - this.stepper = stepper; + public KFitter(final RealVector initialStateEstimate, final RealMatrix initialErrorCovariance, final Propagator propagator, final HashMap materialHashMap) { + this.stateEstimation = initialStateEstimate.copy(); + this.errorCovariance = initialErrorCovariance.copy(); + this.stepper = new Stepper(initialStateEstimate.toArray()); this.propagator = propagator; this.materialHashMap = materialHashMap; } - public void predict(Hit hit, boolean direction) throws Exception { + public void predict(KFHit hit, boolean direction) throws Exception { // Initialization stepper.initialize(direction); Stepper stepper1 = new Stepper(stepper.y); @@ -86,32 +84,12 @@ public void predict(Hit hit, boolean direction) throws Exception { errorCovariance = (transitionMatrix.multiply(errorCovariance.multiply(transitionMatrixT))).add(processNoise); } - public void correct(Hit hit) { - RealVector z; - RealMatrix measurementNoise; - RealMatrix measurementMatrix; - RealVector h; - // check if the hit is the beamline - if (hit.getRadius() < 1) { - // the diagonal elements are the squared errors in r, phi, z - measurementNoise = - new Array2DRowRealMatrix( - new double[][]{ - {vertex_resolutions[0], 0.0000, 0.0000}, - {0.00, 1e10, 0.0000}, - {0.00, 0.0000, vertex_resolutions[1]} - });//3x3 - measurementMatrix = H_beam(stateEstimation);//6x3 - h = h_beam(stateEstimation);//3x1 - 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 - } + public void correct(KFHit hit) { + RealVector z = hit.MeasurementVector(); + RealVector h = hit.ProjectionFunction(stateEstimation); + + RealMatrix measurementNoise = hit.MeasurementNoiseMatrix(); + RealMatrix measurementMatrix = hit.ProjectionMatrix(stateEstimation); RealMatrix measurementMatrixT = measurementMatrix.transpose(); // S = H * P(k) * H' + R @@ -138,6 +116,7 @@ public void correct(Hit hit) { stepper.y = stateEstimation.toArray(); } + // specific to AHDC hits public double residual(Hit hit) { double d = hit.distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); return hit.getDoca()-d; @@ -147,7 +126,7 @@ public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } - private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { + private RealMatrix F(KFHit hit, Stepper stepper1) throws Exception { double[] dfdx = subfunctionF(hit, stepper1, 0); double[] dfdy = subfunctionF(hit, stepper1, 1); @@ -160,7 +139,7 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { {dfdx[0], dfdy[0], dfdz[0], dfdpx[0], dfdpy[0], dfdpz[0]}, {dfdx[1], dfdy[1], dfdz[1], dfdpx[1], dfdpy[1], dfdpz[1]}, {dfdx[2], dfdy[2], dfdz[2], dfdpx[2], dfdpy[2], dfdpz[2]}, {dfdx[3], dfdy[3], dfdz[3], dfdpx[3], dfdpy[3], dfdpz[3]}, {dfdx[4], dfdy[4], dfdz[4], dfdpx[4], dfdpy[4], dfdpz[4]}, {dfdx[5], dfdy[5], dfdz[5], dfdpx[5], dfdpy[5], dfdpz[5]}}); } - double[] subfunctionF(Hit hit, Stepper stepper1, int i) throws Exception { + double[] subfunctionF(KFHit hit, Stepper stepper1, int i) throws Exception { double h = 1e-8;// in mm Stepper stepper_plus = new Stepper(stepper1.y); Stepper stepper_minus = new Stepper(stepper1.y); @@ -184,87 +163,6 @@ private RealMatrix F(Hit hit, Stepper stepper1) throws Exception { return new double[]{dxdi, dydi, dzdi, dpxdi, dpydi, dpzdi}; } - // Measurement matrix in 1x1 dimension: minimize distance - doca - private RealVector h(RealVector x, Hit hit) { - double d = hit.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); - return MatrixUtils.createRealVector(new double[]{d}); - } - - // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) - private RealMatrix H(RealVector x, Hit hit) { - - double ddocadx = subfunctionH(x, hit, 0); - double ddocady = subfunctionH(x, hit, 1); - double ddocadz = subfunctionH(x, hit, 2); - double ddocadpx = subfunctionH(x, hit, 3); - double ddocadpy = subfunctionH(x, hit, 4); - double ddocadpz = subfunctionH(x, hit, 5); - - return MatrixUtils.createRealMatrix(new double[][]{ - {ddocadx, ddocady, ddocadz, ddocadpx, ddocadpy, ddocadpz}}); - } - - double subfunctionH(RealVector x, Hit hit, int i) { - double h = 1e-8;// in mm - RealVector x_plus = x.copy(); - RealVector x_minus = x.copy(); - - x_plus.setEntry(i, x_plus.getEntry(i) + h); - x_minus.setEntry(i, x_minus.getEntry(i) - h); - - double doca_plus = h(x_plus, hit).getEntry(0); - double doca_minus = h(x_minus, hit).getEntry(0); - - return (doca_plus - doca_minus) / (2 * h); - } - - // Measurement matrix for the beamline (Hit_beam) in dimeansion 3x1 - private RealVector h_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - double zz = x.getEntry(2); - - double r = Math.hypot(xx, yy); - double phi = Math.atan2(yy, xx); - - return MatrixUtils.createRealVector(new double[]{r, phi, zz}); - } - - // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) - private RealMatrix H_beam(RealVector x) { - - double xx = x.getEntry(0); - double yy = x.getEntry(1); - - double drdx = (xx) / (Math.hypot(xx, yy)); - double drdy = (yy) / (Math.hypot(xx, yy)); - double drdz = 0.0; - double drdpx = 0.0; - double drdpy = 0.0; - double drdpz = 0.0; - - double dphidx = -(yy) / (xx * xx + yy * yy); - double dphidy = (xx) / (xx * xx + yy * yy); - double dphidz = 0.0; - double dphidpx = 0.0; - double dphidpy = 0.0; - double dphidpz = 0.0; - - double dzdx = 0.0; - double dzdy = 0.0; - double dzdz = 1.0; - double dzdpx = 0.0; - double dzdpy = 0.0; - double dzdpz = 0.0; - - return MatrixUtils.createRealMatrix( - new double[][]{ - {drdx, drdy, drdz, drdpx, drdpy, drdpz}, - {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, - {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} - }); - } public RealVector getStateEstimationVector() { return stateEstimation.copy(); @@ -273,10 +171,8 @@ public RealVector getStateEstimationVector() { public RealMatrix getErrorCovarianceMatrix() { return errorCovariance.copy(); } - - public void setVertexResolution(double[] res) { - vertex_resolutions[0] = res[0]; - vertex_resolutions[1] = res[1]; - } + + /** Return a copy of the stepper. It is like a snapshot of the propagation. */ + public Stepper getStepper() {return new Stepper(stepper);} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java index 268e43862b..0bf0958ff9 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,22 +1,23 @@ package org.jlab.rec.ahdc.KalmanFilter; import java.util.ArrayList; -import java.util.Collections; import java.util.HashMap; +import java.util.Random; import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.RealVector; -import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import org.jlab.clas.tracking.kalmanfilter.Material; -import org.jlab.io.base.DataBank; -import org.jlab.io.base.DataEvent; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.alert.Track.Track; -//import org.apache.commons.math3.linear.RealMatrixFormat; +import org.jlab.detector.calib.utils.DatabaseConstantProvider; +import org.jlab.geom.base.Component; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; +import org.jlab.geom.prim.Point3D; + /** * This is the main routine of the Kalman Filter. The fit is done by a KFitter @@ -31,149 +32,379 @@ public class KalmanFilter { public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;} - public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} + HashMap materialHashMap = MaterialMap.generateMaterials(); private PDGParticle particle; private int Niter = 40; // number of iterations for the Kalman Filter - private boolean IsVtxDefined = false; // implemented but not used yet - private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2} - // mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration - private double clas_alignement = -75; + private double vz_constraint = 0; + private boolean IsVtxDefined = false; + + // mm, they are the misalignement with respect to the AHDC: the are defined in ALERTEngine + private double atof_alignement = 0; - public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { + private int counter = 0; // number of utilisation of the Kalman Filter + + AlertTOFDetector ATOFdet = null; // reference to the ATOF geometry + HashMap> ATOF_hits_predicted = new HashMap<>(); // trackid vs (sector, layer, wedge) + + public void propagation(ArrayList tracks, final double magfield, boolean IsMC) { try { - double vz_constraint = 0; // to be linked to the electron vertex + counter++; // Initialization --------------------------------------------------------------------- final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; - HashMap materialHashMap = MaterialMap.generateMaterials(); - // Recover the vertex of the electron - if (event.hasBank("REC::Particle")) { - DataBank recBank = event.getBank("REC::Particle"); - int row = 0; - while ((!IsVtxDefined) && row < recBank.rows()) { - if (recBank.getInt("pid", row) == 11) { - IsVtxDefined = true; - vz_constraint = 10*recBank.getFloat("vz",row) - (IsMC ? 0 : clas_alignement); // mm - //////////////////////////////////////// - /// compute electron resolution here - /// it depends en p and theta - /// the fine tuning will be done later - /// //////////////////////////////////// - //double px = recBank.getFloat("px",row); - //double py = recBank.getFloat("py",row); - //double pz = recBank.getFloat("pz",row); - //double p = Math.sqrt(px*px+py*py+pz*pz); - //double theta = Math.acos(pz/p); - vertex_resolutions[0] = 0.09; - vertex_resolutions[1] = 64;//4 + 1e10*theta + 1e10*p; - } - row++; - } - } - + + // Loop over tracks for (Track track : tracks) { - // Initialize state vector + /// Initialize state vector double x0 = 0.0; double y0 = 0.0; - double z0 = IsVtxDefined ? vz_constraint : track.get_Z0(); + double z0 = (IsVtxDefined && counter < 2) ? vz_constraint : track.get_Z0(); double px0 = track.get_px(); double py0 = track.get_py(); double pz0 = track.get_pz(); double[] y = new double[]{x0, y0, z0, px0, py0, pz0}; - // Read list of hits + + /// Read list of hits + RadialKFHit beam_hit = track.getBeamlineHit(); ArrayList AHDC_hits = track.getHits(); - Collections.sort(AHDC_hits); // sorted following the compareTo() method in Hit.java + ArrayList ATOF_hits = track.getATOFHits(); - // Start propagation - Stepper stepper = new Stepper(y); + /// Initialize propagator RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B); Propagator propagator = new Propagator(RK4); - // Initialization of the Kalman Fitter - // for the error matrix: first 3 lines in mm^2; last 3 lines in MeV^2 - RealVector initialStateEstimate = new ArrayRealVector(stepper.y); - 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); - if (IsVtxDefined) TrackFitter.setVertexResolution(vertex_resolutions); + /// Initialization of the Kalman Fitter + RealVector initialStateEstimate = new ArrayRealVector(y); + RealMatrix initialErrorCovariance = track.getErrorCovarianceMatrix(); + KFitter TrackFitter = new KFitter(initialStateEstimate, initialErrorCovariance, propagator, materialHashMap); // Loop over number of iterations for (int k = 0; k < Niter; k++) { - // Forward propagation + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { TrackFitter.predict(hit, true); TrackFitter.correct(hit); } - // Backward propagation (last layer to first layer) - for (int i = AHDC_hits.size() - 2; i >= 0; i--) { + + // Take into account ATOF hits + if (!AHDC_hits.isEmpty() && !ATOF_hits.isEmpty()) { + // Forward propagation in ATOF + for (KFHit hit : ATOF_hits) { + TrackFitter.predict(hit, true); + TrackFitter.correct(hit); + } + // Backward propagation in ATOF + for (int i = ATOF_hits.size()-2; i >= 0; i--){ + KFHit hit = ATOF_hits.get(i); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); + } + // Backward propagation to the last AHDC layer + { + Hit hit = AHDC_hits.getLast(); + TrackFitter.predict(hit, false); + TrackFitter.correct(hit); + } + } + + // Backward propagation (from AHDC last layer to first AHDC layer) + for (int i = AHDC_hits.size()-2; i >= 0; i--) { Hit hit = AHDC_hits.get(i); TrackFitter.predict(hit, false); TrackFitter.correct(hit); } - // Backward propagation (first layer to beamline) + + // Backward propagation (from first AHDC layer to beamline) { - Hit hit = new Hit_beam(0, 0, vz_constraint); - TrackFitter.predict(hit, false); - TrackFitter.correct(hit); + TrackFitter.predict(beam_hit, false); + TrackFitter.correct(beam_hit); } - } - - /*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()));*/ + } // end loop over nb. iteration - + /// Output: position and momentum RealVector x_out = TrackFitter.getStateEstimationVector(); track.setPositionAndMomentumVec(x_out.toArray()); + track.setErrorCovarianceMatrix(TrackFitter.getErrorCovarianceMatrix()); - // Post fit propagation (no correction) to set the residuals - KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap); + /// Post fit propagation (no correction) + KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Propagator(RK4), materialHashMap); + + // Forward propagation in AHDC for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); - if( hit.getId()>0){ // for the beamline the hit id is 0, so we only look at AHDC hits - hit.setResidual(PostFitPropagator.residual(hit)); - } + hit.setResidual(PostFitPropagator.residual(hit)); // output: residual } - - // Fill track and hit bank + + // Fill values for AHDC::kftrack // TO DO : s and p_drift have to be checked to be sure they represent what we want - double s = PostFitPropagator.stepper.sTot; - double p_drift = PostFitPropagator.stepper.p(); - int sum_adc = 0; + Stepper current_stepper = PostFitPropagator.getStepper(); + double s = current_stepper.sTot; + double p_drift = current_stepper.p(); double sum_residuals = 0; double chi2 = 0; for (Hit hit : AHDC_hits) { - sum_adc += hit.getADC(); sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0); + chi2 += Math.pow(hit.getResidual(),2)/hit.MeasurementNoiseMatrix().getEntry(0,0); } - track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); track.set_chi2(chi2/(AHDC_hits.size()-3)); track.set_p_drift(p_drift); - track.set_dEdx(sum_adc/s); track.set_path(s); - track.set_n_hits(AHDC_hits.size()); + + /// At the end of the PostFit propagation towards the last layer of the AHDC + /// we project the track on the ATOF surface without any AI ATOF hit indication + if (ATOFdet != null) { + // Projection towards the ATOF surfaces + // R1 : radius of the lower surface of an ATOF bar + // R2 : radius of the upper surface of an ATOF bar = lower surface of an ATOF wedge + // R3 : radius of the upper surface of an ATOF wedge + // Illustration of an ATOF component : p0 and p3 are the points on the lower surface; p1 and p2 are the points on the upper surface + // p2 ---------- p1 + // \ / + // \ / + // p3 ------ p0 + Point3D pR1 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(0); + Point3D pR2 = ATOFdet.getSector(0).getSuperlayer(0).getLayer(0).getComponent(10).getVolumePoint(2); + Point3D pR3 = ATOFdet.getSector(0).getSuperlayer(1).getLayer(0).getComponent(0).getVolumePoint(2); + double R1 = Math.hypot(pR1.x(), pR1.y()); + double R2 = Math.hypot(pR2.x(), pR2.y()); + double R3 = Math.hypot(pR3.x(), pR3.y()); + + /// From last AHDC hit to surface R1 + RadialSurfaceKFHit hitR1 = new RadialSurfaceKFHit(R1); + PostFitPropagator.predict(hitR1, true); + Stepper stepperR1 = PostFitPropagator.getStepper(); + track.set_ATOF_S1_stepper(stepperR1); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR1.r() - R1) < 1.5*stepperR1.h) { + track.set_ATOF_region(1); // indicate if we reach this surface or not + } + + /// From surface R1 to surface R2 + RadialSurfaceKFHit hitR2 = new RadialSurfaceKFHit(R2); + PostFitPropagator.predict(hitR2, true); + Stepper stepperR2 = PostFitPropagator.getStepper(); + track.set_ATOF_S2_stepper(stepperR2); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR2.r() - R2) < 1.5*stepperR2.h) { + track.set_ATOF_region(2); // indicate if we reach this surface or not + } + // predict the wedge + double[] vecR2 = PostFitPropagator.getStateEstimationVector().toArray(); + Point3D posR2 = new Point3D(vecR2[0], vecR2[1], vecR2[2]); + posR2.translateXYZ(0, 0, atof_alignement); + ATOF_hits_predicted.put(track.get_trackId(), predict_adjacent_wedges(ATOFdet, posR2)); + + + /// From surface R2 to surface R3 + RadialSurfaceKFHit hitR3 = new RadialSurfaceKFHit(R3); + PostFitPropagator.predict(hitR3, true); + Stepper stepperR3 = PostFitPropagator.getStepper(); + track.set_ATOF_S3_stepper(stepperR3); // to store : x,y,z, path, p in AHDC::kftrack + if (Math.abs(stepperR3.r() - R3) < 1.5*stepperR3.h) { + track.set_ATOF_region(3); // indicate if we reach this surface or not + } + + } // end propagation towards ATOF surface + }//end of loop on track candidates } catch (Exception e) { - //e.printStackTrace(); - //System.out.println("======> Kalman Filter Error"); + e.printStackTrace(); + //System.out.println("Error in Kalman Filter"); + } + } + public void set_Niter(int Niter) {this.Niter = Niter;} + public int get_Niter() {return this.Niter;} + public void set_particle(PDGParticle particle) {this.particle = particle;} + public PDGParticle get_particle() {return this.particle;} + public HashMap> get_ATOF_hits_predicted() {return this.ATOF_hits_predicted;} + + // Test + public static void main(String[] args) { + // ATOF detector + AlertTOFDetector atof = (new AlertTOFFactory()).createDetectorCLAS(new DatabaseConstantProvider()); + + int Npts = 1000; + int counter = 0; + int err = 0; + for (int i = 0; i < Npts; i++) { + Random rand = new Random(); + int true_sector = rand.nextInt(15); + int true_layer = rand.nextInt(4); + int true_wedge = rand.nextInt(10); + + Component comp = atof.getSector(true_sector).getSuperlayer(1).getLayer(true_layer).getComponent(true_wedge); + // top face + Point3D p0 = comp.getVolumePoint(0); + Point3D p1 = comp.getVolumePoint(1); + Point3D p2 = comp.getVolumePoint(2); + Point3D p3 = comp.getVolumePoint(3); + // bottom face + Point3D p4 = comp.getVolumePoint(4); + Point3D p5 = comp.getVolumePoint(5); + Point3D p6 = comp.getVolumePoint(6); + Point3D p7 = comp.getVolumePoint(7); + + // Random point int he current wedge volume + //Point3D pt = p0.lerp(p7, Math.random()); + //Point3D pt = comp.getMidpoint(); + double t0 = rand.nextDouble(1); + double t1 = rand.nextDouble(1-t0); + double t2 = rand.nextDouble(1-t0-t1); + double t3 = 1-t0-t1-t2; + //System.out.printf("t0 + t1 + t2 + t3 = %f\n", t0+t1+t2+t3); + double t4 = rand.nextDouble(1); + double t5 = rand.nextDouble(1-t4); + double t6 = rand.nextDouble(1-t4-t5); + double t7 = 1-t4-t5-t6; + //System.out.printf("t4 + t5 + t6 + t7 = %f\n", t4+t5+t6+t7); + double x_top = t0*p0.x() + t1*p1.x() + t2*p2.x() + t3*p3.x(); + double y_top = t0*p0.y() + t1*p1.y() + t2*p2.y() + t3*p3.y(); + double x_bot = t4*p4.x() + t5*p5.x() + t6*p6.x() + t7*p7.x(); + double y_bot = t4*p4.y() + t5*p5.y() + t6*p6.y() + t7*p7.y(); + Point3D pt_top = new Point3D(x_top, y_top, p0.z()); + Point3D pt_bot = new Point3D(x_bot, y_bot, p4.z()); + Point3D pt = pt_top.lerp(pt_bot, Math.random()); + //System.out.printf("distance from midpoint : %f\n", comp.getMidpoint().distance(pt)); + + // Test the algoritm + int[] res = KalmanFilter.predict_wedge(atof, pt); + + if (res[0] == true_sector && res[1] == true_layer && res[2] == true_wedge) { + counter++; + } else { + err++; + System.out.printf("%d) Initial wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, true_sector, true_layer, true_wedge); + System.out.printf("%d) Predicted wedge : sector (% 2d) layer (% 2d) wedge (% 2d)\n", err, res[0], res[1], res[2]); + } + } + System.out.printf("Nb of testing : %d\n", Npts); + System.out.printf("Nb of success : %d (%.2f %%)\n", counter, 100.0*counter/Npts); + } + /** + * @param pt is defined in the center of the ATOF + */ + static public int[] predict_wedge(AlertTOFDetector atof, Point3D pt) { + // find the wedge + int wedge = -1; + double dz = 1e10; + for (int c = 0; c < atof.getSector(0).getSuperlayer(1).getLayer(0).getNumComponents(); c++) { + Point3D midpoint = atof.getSector(0).getSuperlayer(1).getLayer(0).getComponent(c).getMidpoint(); + if (Math.abs(midpoint.z()-pt.z()) < dz) { + dz = Math.abs(midpoint.z()-pt.z()); + wedge = c; + } + } + if (wedge == -1) return null; + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(1).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(1).getLayer(l).getComponent(wedge).getMidpoint(); + if (midpoint.distance(pt) < d) { + d = midpoint.distance(pt); + sector = s; + layer = l; + } + } + } + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, wedge}; } } - void set_Niter(int Niter) {this.Niter = Niter;} - void set_particle(PDGParticle particle) {this.particle = particle;} + + /** + * @param pt is defined in the center of the ATOF + */ + static public int[] predict_bar(AlertTOFDetector atof, Point3D pt) { + // find sector and layer + int sector = -1; + int layer = -1; + double d = 1e10; + for (int s = 0; s < atof.getNumSectors(); s++) { + for (int l = 0; l < atof.getSector(s).getSuperlayer(0).getNumLayers(); l++) { + Point3D midpoint = atof.getSector(s).getSuperlayer(0).getLayer(l).getComponent(10).getMidpoint(); + double distance = midpoint.vectorTo(pt).rho(); + if (distance < d) { + d = distance; + sector = s; + layer = l; + } + } + } + if (sector == -1 || layer == -1) { + return null; + } else { + return new int[] {sector, layer, 10}; + } + } + + public ArrayList predict_adjacent_wedges(AlertTOFDetector atof, Point3D pt) { + int[] closest_wedge_id = predict_wedge(atof, pt); + if (closest_wedge_id == null) { + return new ArrayList(); // an empty list + } + + int sector = closest_wedge_id[0]; + int layer = closest_wedge_id[1]; + int wedge = closest_wedge_id[2]; + + // find adjacent layer and sector + int sector_plus = sector; + int sector_minus = sector; + int layer_plus = layer+1; + int layer_minus = layer-1; + if (layer == 0) { + sector_plus = sector; + sector_minus = Math.floorMod(sector-1, 15); + layer_plus = layer+1; + layer_minus = 3; + } + else if (layer == 3) { + sector_plus = Math.floorMod(sector+1, 15); + sector_minus = sector; + layer_plus = 0; + layer_minus = layer-1; + } + // Here are all the adjacents wedges (maximum 9) + ArrayList listOfWedges = new ArrayList<>(); + listOfWedges.add(new int[]{sector, layer, wedge}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge}); + if (wedge-1 >= 0) { + listOfWedges.add(new int[]{sector, layer, wedge-1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge-1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge-1}); + } + if (wedge+1 <= 9) { + listOfWedges.add(new int[]{sector, layer, wedge+1}); + listOfWedges.add(new int[]{sector_plus, layer_plus, wedge+1}); + listOfWedges.add(new int[]{sector_minus, layer_minus, wedge+1}); + } + + // order the list following the distance to the initial point + listOfWedges.sort((a,b) -> { + double da = atof.getSector(a[0]).getSuperlayer(1).getLayer(a[1]).getComponent(a[2]).getMidpoint().distance(pt); + double db = atof.getSector(b[0]).getSuperlayer(1).getLayer(b[1]).getComponent(b[2]).getMidpoint().distance(pt); + return Double.compare(da, db); + }); + + return listOfWedges; + } + + public void set_ATOF_detector(AlertTOFDetector atof) { this.ATOFdet = atof;} + public void set_atof_alignement(double _shift) {this.atof_alignement = _shift;} + public void set_vz_constraint(double _vz) {this.vz_constraint = _vz;} + public void set_vertex_flag(boolean _flag) {this.IsVtxDefined = _flag;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java index e85021addb..6f4f58a835 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/Propagator.java @@ -1,13 +1,10 @@ package org.jlab.rec.ahdc.KalmanFilter; -//import java.util.Arrays; - import org.apache.commons.math3.linear.ArrayRealVector; import org.apache.commons.math3.linear.RealVector; import org.jlab.geom.prim.Point3D; import org.jlab.clas.tracking.kalmanfilter.Material; import java.util.HashMap; -import org.jlab.rec.ahdc.Hit.Hit; /** * Is responsible of the propagation @@ -27,7 +24,7 @@ public Propagator(RungeKutta4 rungeKutta4) { } // Propagate the stepper toward the next hit - void propagate(Stepper stepper, Hit hit, HashMap materialHashMap) { + void propagate(Stepper stepper, KFHit hit, HashMap materialHashMap) { // Do not allow more than 10000 steps (very critical cases) final int maxNbOfStep = 10000; @@ -95,9 +92,8 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } else { // the distance between the step is not so big but the distance with respect to the hit starts to increase - // so, go back to the previous step (the best we have) and stop the propagation, set the default step size + // so, go back to the previous step (the best we have) and stop the propagation stepper.copyContent(prev_stepper); - stepper.h = 0.5; break; } } @@ -161,7 +157,7 @@ void propagate(Stepper stepper, Hit hit, HashMap materialHashM } - public RealVector f(Stepper stepper, Hit hit, HashMap materialHashMap) { + public RealVector f(Stepper stepper, KFHit hit, HashMap materialHashMap) { propagate(stepper, hit, materialHashMap); return new ArrayRealVector(stepper.y); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java new file mode 100644 index 0000000000..220118255b --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -0,0 +1,105 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; + +/** + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ +public class RadialKFHit implements KFHit { + + private double r,phi,z; + Line3D line; + RealMatrix measurementNoise = null; + + public RadialKFHit(double x, double y , double z) { + this.z = z; + this.r = Math.hypot(x,y); + this.phi = Math.atan2(y,x); + line = new Line3D(x,y,0,x,y,1); // a line parallel to the beam axis + } + + @Override + public RealVector MeasurementVector() { + return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); + } + + @Override + public RealMatrix MeasurementNoiseMatrix() { + return measurementNoise; + } + + public void setMeasurementNoise(RealMatrix measurementNoise) { + this.measurementNoise= measurementNoise; + } + + @Override + public double distance(Point3D point3D) { + return this.line.distance(point3D).length(); + } + + @Override + public double getRadius() { + return r; + } + + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + double zz = x.getEntry(2); + + double r = Math.hypot(xx, yy); + double phi = Math.atan2(yy, xx); + + return MatrixUtils.createRealVector(new double[]{r, phi, zz}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + double dphidx = -(yy) / (xx * xx + yy * yy); + double dphidy = (xx) / (xx * xx + yy * yy); + double dphidz = 0.0; + double dphidpx = 0.0; + double dphidpy = 0.0; + double dphidpz = 0.0; + + double dzdx = 0.0; + double dzdy = 0.0; + double dzdz = 1.0; + double dzdpx = 0.0; + double dzdpy = 0.0; + double dzdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz}, + {dphidx, dphidy, dphidz, dphidpx, dphidpy, dphidpz}, + {dzdx, dzdy, dzdz, dzdpx, dzdpy, dzdpz} + }); + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java new file mode 100644 index 0000000000..6e633877bb --- /dev/null +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java @@ -0,0 +1,79 @@ +package org.jlab.rec.ahdc.KalmanFilter; + +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RealVector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; + +/** + * Implement a hit for which the state vector is a 3x1 matrix (r, phi, z) + * e.g beamline, ATOF wedges/bars + * + * @author Mathieu Ouillon + * @author Éric Fuchey + * @author Felix Touchte Codjo + */ +public class RadialSurfaceKFHit implements KFHit { + + private double r; + + public RadialSurfaceKFHit(double r) { + this.r = r; + } + + @Override + public RealVector MeasurementVector() { + return new ArrayRealVector(new double[] {this.r}); + } + + @Override + public RealMatrix MeasurementNoiseMatrix() { + return new Array2DRowRealMatrix(new double[][]{{1e-8}}); + } + + @Override + public double distance(Point3D point3D) { + return Math.abs(this.r - Math.hypot(point3D.x(), point3D.y())); + } + + @Override + public double getRadius() { + return r; + } + + // Projection function + @Override + public RealVector ProjectionFunction(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double r = Math.hypot(xx, yy); + + return MatrixUtils.createRealVector(new double[]{r}); + } + + // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) + @Override + public RealMatrix ProjectionMatrix(RealVector x) { + + double xx = x.getEntry(0); + double yy = x.getEntry(1); + + double drdx = (xx) / (Math.hypot(xx, yy)); + double drdy = (yy) / (Math.hypot(xx, yy)); + double drdz = 0.0; + double drdpx = 0.0; + double drdpy = 0.0; + double drdpz = 0.0; + + return MatrixUtils.createRealMatrix( + new double[][]{ + {drdx, drdy, drdz, drdpx, drdpy, drdpz} + }); + } + +} diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java index e108ee3493..09cfb97b81 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/alert/Track/Track.java @@ -1,9 +1,13 @@ package org.jlab.rec.alert.Track; import org.jlab.rec.alert.AI.InterCluster; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.rec.ahdc.AHDCCluster.AHDCCluster; import org.jlab.rec.ahdc.HelixFit.HelixFitObject; import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; +import org.jlab.rec.ahdc.KalmanFilter.Stepper; import java.util.ArrayList; import java.util.List; @@ -140,4 +144,67 @@ public double get_dEdx() { public int get_predicted_ATOF_sector() { return candidate.get_predicted_ATOF_sector(); } public int get_predicted_ATOF_layer() { return candidate.get_predicted_ATOF_layer(); } public int get_predicted_ATOF_wedge() { return candidate.get_predicted_ATOF_wedge(); } + + // Recent extension of the Kalman Filter + /** + *

Error covariance matrix of the measurment around the beamline

+ * + *

This quite useful when one plan to run the Kalman filter several times. This object stores the current status of the error covariance matrix at the end of the backward propagation of an iteration.

+ * + *

First 3 lines in mm^2; last 3 lines in MeV^2 (in the beamline)

+ */ + RealMatrix errorCovarianceMatrix = MatrixUtils.createRealMatrix(new double[][]{ + {50 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 50 , 0.0 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 900 , 0.0 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 100 , 0.0 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 100 , 0.0}, + {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 900}}); + + /** + * Get the starting error covariance matrix {@link #errorCovarianceMatrix} (around the beamline) for the current KF iteration. + */ + public RealMatrix getErrorCovarianceMatrix() {return errorCovarianceMatrix;} + + private ArrayList ATOF_hits = new ArrayList<>(); + private RadialKFHit beamline_hit = null; + + public ArrayList getATOFHits() { + return this.ATOF_hits; + } + + public RadialKFHit getBeamlineHit() { + return this.beamline_hit; + } + + public void setATOFHits(ArrayList _ATOF_hits) {this.ATOF_hits = _ATOF_hits;} + + public void setBeamlineHit(RadialKFHit _beamline_hit) {this.beamline_hit = _beamline_hit;} + + /** + * Update the error covariance matrix {@link #errorCovarianceMatrix} + */ + public void setErrorCovarianceMatrix(RealMatrix errorCovarianceMatrix) {this.errorCovarianceMatrix = errorCovarianceMatrix;} + + // Position and momentum when the track crosses the ATOF surface + // S1 : lower surface of an ATOF bar + // S2 : upper surface of an ATOF bar = lower surface of an ATOF wedge + // S3 : upper surface of an ATOF wedge + Stepper ATOF_S1_stepper; + Stepper ATOF_S2_stepper; + Stepper ATOF_S3_stepper; + double ATOF_S1_radius; + double ATOF_S2_radius; + double ATOF_S3_radius; + int ATOF_region = 0; // is n if the trach reaches Sn, 0 otherwise (i.e does not reach S1) + + // Projection of the Track on the ATOF surfaces + public void set_ATOF_S1_stepper(Stepper _stepper) {this.ATOF_S1_stepper = _stepper;} + public void set_ATOF_S2_stepper(Stepper _stepper) {this.ATOF_S2_stepper = _stepper;} + public void set_ATOF_S3_stepper(Stepper _stepper) {this.ATOF_S3_stepper = _stepper;} + public void set_ATOF_region(int _n) {this.ATOF_region = _n;} + public Stepper get_ATOF_S1_stepper() {return this.ATOF_S1_stepper;} + public Stepper get_ATOF_S2_stepper() {return this.ATOF_S2_stepper;} + public Stepper get_ATOF_S3_stepper() {return this.ATOF_S3_stepper;} + public int get_ATOF_region() {return this.ATOF_region;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 9b2ba0f773..5167e70774 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -8,6 +8,7 @@ import java.util.HashSet; import java.util.Set; import java.util.concurrent.atomic.AtomicInteger; +import java.util.Iterator; import java.util.HashMap; import java.util.Map; @@ -17,6 +18,7 @@ import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.utils.groups.IndexedTable; import org.jlab.geom.base.Detector; +import org.jlab.geom.detector.alert.ATOF.AlertTOFDetector; import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; import org.jlab.io.base.DataBank; import org.jlab.io.base.DataEvent; @@ -33,6 +35,7 @@ import org.jlab.rec.ahdc.DocaCluster.DocaClusterRefiner; import org.jlab.rec.ahdc.HelixFit.HelixFitJava; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; +import org.jlab.rec.ahdc.KalmanFilter.RadialKFHit; import org.jlab.rec.ahdc.TrackFindingMode; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.PreCluster.PreClusterFinder; @@ -48,6 +51,8 @@ import org.jlab.rec.alert.Track.AtofHitStub; import org.jlab.rec.alert.Track.Track; import org.jlab.rec.alert.Track.TrackCandidate; +import org.apache.commons.math3.linear.Array2DRowRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; import org.jlab.clas.pdg.PDGDatabase; import org.jlab.clas.pdg.PDGParticle; import java.util.List; @@ -99,7 +104,7 @@ public class ALERTEngine extends ReconstructionEngine { */ private RecoBankWriter rbc; static final Logger LOGGER = Logger.getLogger(ALERTEngine.class.getName()); - Detector ATOF; // ALERT ATOF detector + private AlertTOFDetector ATOF; // ALERT ATOF detector private AlertDCDetector AHDC; // ALERT AHDC detector /** @@ -521,8 +526,42 @@ public boolean processDataEventUser(DataEvent event) { /// Kalmam Filter /// /////////////////////////////////////// - // read the list of tracks/hits from the banks AHDC::track and AHDC::hits + /// Pre conditions if (!event.hasBank("AHDC::track")) {return false;} + if (!event.hasBank("AHDC::hits")) {return false;} + + /// tmp: misalignement with respect to the center of the AHDC (mm) + double clas_alignement = +75; + double atof_alignement = 0; + + /// Read the electron vertex + double vz_electron = 0; + double[] vz_error2 = {0.09, 1e10, 1e10}; // mm^2, radians^2, mm^2, error on r, phi, z + boolean IsVertexDefined = false; + if (event.hasBank("REC::Particle")) { + DataBank recBank = event.getBank("REC::Particle"); + for (int row = 0; row < recBank.rows(); row++) { + if (recBank.getInt("pid", row) == 11) { + vz_electron = 10*recBank.getFloat("vz",row); // conversion in mm + IsVertexDefined = true; + + + //double px = recBank.getFloat("px",row); + //double py = recBank.getFloat("py",row); + //double pz = recBank.getFloat("pz",row); + //double p = Math.sqrt(px*px+py*py+pz*pz); + //double theta = Math.acos(pz/p); + + // set the resolutions on r and z! to be done + vz_error2[0] = 0.09; // should depend on p and theta + vz_error2[2] = 64; // should depend on p and theta + + break; // only look at the first electron + } + } + } + + /// Read the list of tracks/hits from the banks AHDC::track and AHDC::hits DataBank trackBank = event.getBank("AHDC::track"); DataBank hitBank = event.getBank("AHDC::hits"); ArrayList AHDC_tracks = new ArrayList<>(); @@ -568,25 +607,156 @@ public boolean processDataEventUser(DataEvent event) { newTrack.set_trackId(trackid); AHDC_tracks.add(newTrack); } - // intialise the Kalman Filter + + /// Associate the electron vertex (the beamline hit) to each track + boolean IsMC = event.hasBank("MC::Particle"); + double vz_constraint = vz_electron + (IsMC ? 0 : clas_alignement); // we don't have the misalignment in simulation + for (Track track : AHDC_tracks) { + RadialKFHit hit_beam = new RadialKFHit(0, 0, vz_constraint); + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {vz_error2[0], 0.0000 , 0.0000}, + {0.0000 , vz_error2[1], 0.0000}, + {0.0000 , 0.0000 , vz_error2[2]} + });//3x3; + hit_beam.setMeasurementNoise(measurementNoise); + track.setBeamlineHit(hit_beam); + } + + /// Look for ATOF wedge hits predicted by the AI + HashMap map_ATOF_hits = new HashMap<>(); + for (Pair pair : matched_ATOF_hit_id) { + int trackid = pair.getKey(); + int atofid = pair.getValue(); + if (trackid > 0 && atofid > 0) { + // recover the wedge + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getShort("id", row) == atofid) { + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + z += atof_alignement; // there is a shift between AHDC and ATOF (still don't know why) ! + RadialKFHit hit = new RadialKFHit(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.setMeasurementNoise(measurementNoise); + map_ATOF_hits.put(trackid, hit); + } + } + } + } + + /// Associate the ATOF hits to each track + for (Track track : AHDC_tracks) { + RadialKFHit hit = map_ATOF_hits.get(track.get_trackId()); + ArrayList list = new ArrayList<>(); + if (hit != null) list.add(hit); // for now, we only consider one hit in the ATOF + track.setATOFHits(list); + } + + /// Intialise the Kalman Filter double magfieldfactor = runBank.getFloat("solenoid", 0); double magfield = 50*magfieldfactor; - boolean IsMC = event.hasBank("MC::Particle"); PDGParticle proton = PDGDatabase.getParticleById(2212); - int Niter = 40; + int Niter = 25; KalmanFilter KF = new KalmanFilter(proton, Niter); - /////////////////////////////////////////////////////// - // first propagation : each AHDC_tracks will be fitted - /////////////////////////////////////////////////////// - KF.propagation(AHDC_tracks, event, magfield, IsMC); - ///////////////////////////////////////////// - // write the AHDC::kftrack bank in the event - ///////////////////////////////////////////// + KF.set_ATOF_detector(null); + //KF.set_ATOF_detector(ATOF); // Reference the ATOF geometry in the Kalman Filter + KF.set_atof_alignement(atof_alignement); + KF.set_vz_constraint(vz_constraint); + KF.set_vertex_flag(IsVertexDefined); + + /// Do a first propagation + KF.propagation(AHDC_tracks, magfield, IsMC); + + /// Look at the new ATOF hits predicted after projection of the track on the lower surface of the ATOF wedges + HashMap> ATOF_hits_predicted = KF.get_ATOF_hits_predicted(); + for (Track track : AHDC_tracks) { + int trackid = track.get_trackId(); + ArrayList possible_wedges = ATOF_hits_predicted.get(trackid); + boolean IsHitSelected = false; + if (possible_wedges != null) { + for (int[] id : possible_wedges) { + int sector = id[0]; + int layer = id[1]; + int wedge = id[2]; + // check if this hit exist in ATOF::hits + for (int row = 0; row < bank_ATOFHits.rows(); row++) { + if (bank_ATOFHits.getInt("sector", row) == sector && bank_ATOFHits.getInt("layer", row) == layer && bank_ATOFHits.getInt("component", row) == wedge) { + // create a RadialKFHit + double x = bank_ATOFHits.getFloat("x", row); + double y = bank_ATOFHits.getFloat("y", row); + double z = bank_ATOFHits.getFloat("z", row); + z += atof_alignement; // there is a shift between AHDC and ATOF (still don't know why) ! + RadialKFHit hit = new RadialKFHit(x, y, z); + // error on r + double wedge_width = 20; //mm + double dr2 = Math.pow(wedge_width, 2)/12; // mm^2 + // error on phi + double open_angle = Math.toRadians(6); // deg + double dphi2 = Math.pow(open_angle, 2)/12; + // error on z + double wedge_length = 27.7; //mm + double dz2 = Math.pow(wedge_length, 2)/12; + + RealMatrix measurementNoise = new Array2DRowRealMatrix( + new double[][]{ + {dr2, 0.0000, 0.0000}, + {0.00, dphi2, 0.0000}, + {0.00, 0.0000, dz2} + });//3x3; + hit.setMeasurementNoise(measurementNoise); + + ArrayList list = new ArrayList<>(); + list.add(hit); // for now, we only consider one hit in the ATOF + track.setATOFHits(list); // update the list of the ATOF hit (i.e override AI ATOF hit if it exist) + + IsHitSelected = true; + break; + } // end id matching + } // end loop over atof hit + if (IsHitSelected) break; + } + } + } + + /// Clean AHDC bad hits + double sigma = 0.5; // mm + for (Track track : AHDC_tracks) { + ArrayList AHDC_hits = track.getHits(); + Iterator it = AHDC_hits.iterator(); + while (it.hasNext()) { + Hit hit = it.next(); + if (Math.abs(hit.getResidual()) > 3*sigma) { + it.remove(); + } + } + } + + /// Second propagation : each AHDC_tracks will be fitted + KF.set_Niter(15); + KF.propagation(AHDC_tracks, magfield, IsMC); + + /// Write the AHDC::kftrack bank in the event event.removeBank("AHDC::kftrack"); org.jlab.rec.ahdc.Banks.RecoBankWriter ahdc_writer = new org.jlab.rec.ahdc.Banks.RecoBankWriter(); DataBank recoKFTracksBank = ahdc_writer.fillAHDCKFTrackBank(event, AHDC_tracks); event.appendBank(recoKFTracksBank); - // update the AHDC::hits bank : fill the residuals + /// Update the AHDC::hits bank : fill the residuals event.removeBank("AHDC::hits"); ArrayList AHDC_hits = new ArrayList<>(); for (Track track : AHDC_tracks) { From aae037e7ac10ffe9b408def37ba5259ecb4e2bd4 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Tue, 9 Jun 2026 13:40:40 -0400 Subject: [PATCH 4/8] activate atof hit --- .../src/main/java/org/jlab/service/alert/ALERTEngine.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index 5167e70774..7e9428ff74 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -674,8 +674,8 @@ public boolean processDataEventUser(DataEvent event) { PDGParticle proton = PDGDatabase.getParticleById(2212); int Niter = 25; KalmanFilter KF = new KalmanFilter(proton, Niter); - KF.set_ATOF_detector(null); - //KF.set_ATOF_detector(ATOF); // Reference the ATOF geometry in the Kalman Filter + //KF.set_ATOF_detector(null); + KF.set_ATOF_detector(ATOF); // Reference the ATOF geometry in the Kalman Filter KF.set_atof_alignement(atof_alignement); KF.set_vz_constraint(vz_constraint); KF.set_vertex_flag(IsVertexDefined); From 555c315f91aa1b16498d434c42b6db821ffd29e9 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 10 Jun 2026 05:07:16 -0400 Subject: [PATCH 5/8] rename methods --- .../src/main/java/org/jlab/rec/ahdc/Hit/Hit.java | 12 ++++++------ .../java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java | 10 +++++----- .../java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 8 ++++---- .../org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 2 +- .../org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java | 8 ++++---- .../rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java | 8 ++++---- 6 files changed, 24 insertions(+), 24 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 18dd1bcad5..cd97e04967 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 @@ -193,12 +193,12 @@ public void setTrackId(int _trackId) { } @Override - public RealVector MeasurementVector() { + public RealVector getMeasurementVector() { return new ArrayRealVector(new double[]{this.doca}); } @Override - public RealMatrix MeasurementNoiseMatrix() { + public RealMatrix getMeasurementNoiseMatrix() { double mean_error = 0.471; // mm (no difference between adc and time) double error_on_adc = (1.15146*raw_adc + 437.63)/(3.21187*raw_adc + 878.855); // mm double error_on_time = (0.4423*time + 13.7215)/(0.846038*time + 31.9867); // mm @@ -210,14 +210,14 @@ public RealMatrix MeasurementNoiseMatrix() { // Projection function @Override - public RealVector ProjectionFunction(RealVector x) { + public RealVector getProjectionFunction(RealVector x) { double d = this.distance(new Point3D(x.getEntry(0), x.getEntry(1), x.getEntry(2))); return MatrixUtils.createRealVector(new double[]{d}); } // Jacobian matrix of the measurement with respect to (x, y, z, px, py, pz) @Override - public RealMatrix ProjectionMatrix(RealVector x) { + public RealMatrix getProjectionMatrix(RealVector x) { double ddocadx = partialProjectionMatrix(x, 0); double ddocady = partialProjectionMatrix(x, 1); @@ -238,8 +238,8 @@ private double partialProjectionMatrix(RealVector x, int i) { x_plus.setEntry(i, x_plus.getEntry(i) + h); x_minus.setEntry(i, x_minus.getEntry(i) - h); - double doca_plus = this.ProjectionFunction(x_plus).getEntry(0); - double doca_minus = this.ProjectionFunction(x_minus).getEntry(0); + double doca_plus = this.getProjectionFunction(x_plus).getEntry(0); + double doca_minus = this.getProjectionFunction(x_minus).getEntry(0); return (doca_plus - doca_minus) / (2 * h); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java index 3a2e693192..399bd6f71c 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFHit.java @@ -12,11 +12,11 @@ public interface KFHit { public double distance(Point3D point3D); public double getRadius(); /** Return the measurement encoded in this KFHit */ - public RealVector MeasurementVector(); + public RealVector getMeasurementVector(); /** Return the measurement noise matrix for this this KFHit */ - public RealMatrix MeasurementNoiseMatrix(); + public RealMatrix getMeasurementNoiseMatrix(); /** Compute the measure for a given state vector */ - public RealVector ProjectionFunction(RealVector x); - /** Compute the Jacobian matrix of the {@link #ProjectionFunction(RealVector)} with respect of the components of the state vector */ - public RealMatrix ProjectionMatrix(RealVector x); + public RealVector getProjectionFunction(RealVector x); + /** Compute the Jacobian matrix of the {@link #getProjectionFunction(RealVector)} with respect of the components of the state vector */ + public RealMatrix getProjectionMatrix(RealVector x); } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KFitter.java index 8fb0ecaafd..12eee314bd 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 @@ -85,11 +85,11 @@ public void predict(KFHit hit, boolean direction) throws Exception { } public void correct(KFHit hit) { - RealVector z = hit.MeasurementVector(); - RealVector h = hit.ProjectionFunction(stateEstimation); + RealVector z = hit.getMeasurementVector(); + RealVector h = hit.getProjectionFunction(stateEstimation); - RealMatrix measurementNoise = hit.MeasurementNoiseMatrix(); - RealMatrix measurementMatrix = hit.ProjectionMatrix(stateEstimation); + RealMatrix measurementNoise = hit.getMeasurementNoiseMatrix(); + RealMatrix measurementMatrix = hit.getProjectionMatrix(stateEstimation); RealMatrix measurementMatrixT = measurementMatrix.transpose(); // S = H * P(k) * H' + R 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 0bf0958ff9..371908ec12 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 @@ -153,7 +153,7 @@ public void propagation(ArrayList tracks, final double magfield, boolean double chi2 = 0; for (Hit hit : AHDC_hits) { sum_residuals += hit.getResidual(); - chi2 += Math.pow(hit.getResidual(),2)/hit.MeasurementNoiseMatrix().getEntry(0,0); + chi2 += Math.pow(hit.getResidual(),2)/hit.getMeasurementNoiseMatrix().getEntry(0,0); } track.set_sum_residuals(sum_residuals); track.set_chi2(chi2/(AHDC_hits.size()-3)); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java index 220118255b..3f6e7c60f3 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialKFHit.java @@ -29,12 +29,12 @@ public RadialKFHit(double x, double y , double z) { } @Override - public RealVector MeasurementVector() { + public RealVector getMeasurementVector() { return new ArrayRealVector(new double[] {this.r, this.phi, this.z}); } @Override - public RealMatrix MeasurementNoiseMatrix() { + public RealMatrix getMeasurementNoiseMatrix() { return measurementNoise; } @@ -54,7 +54,7 @@ public double getRadius() { // Projection function @Override - public RealVector ProjectionFunction(RealVector x) { + public RealVector getProjectionFunction(RealVector x) { double xx = x.getEntry(0); double yy = x.getEntry(1); @@ -68,7 +68,7 @@ public RealVector ProjectionFunction(RealVector x) { // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) @Override - public RealMatrix ProjectionMatrix(RealVector x) { + public RealMatrix getProjectionMatrix(RealVector x) { double xx = x.getEntry(0); double yy = x.getEntry(1); diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java index 6e633877bb..736f46cba6 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/RadialSurfaceKFHit.java @@ -25,12 +25,12 @@ public RadialSurfaceKFHit(double r) { } @Override - public RealVector MeasurementVector() { + public RealVector getMeasurementVector() { return new ArrayRealVector(new double[] {this.r}); } @Override - public RealMatrix MeasurementNoiseMatrix() { + public RealMatrix getMeasurementNoiseMatrix() { return new Array2DRowRealMatrix(new double[][]{{1e-8}}); } @@ -46,7 +46,7 @@ public double getRadius() { // Projection function @Override - public RealVector ProjectionFunction(RealVector x) { + public RealVector getProjectionFunction(RealVector x) { double xx = x.getEntry(0); double yy = x.getEntry(1); @@ -58,7 +58,7 @@ public RealVector ProjectionFunction(RealVector x) { // Jacobian matrix of the measurement for the beamline with respect to (x, y, z, px, py, pz) @Override - public RealMatrix ProjectionMatrix(RealVector x) { + public RealMatrix getProjectionMatrix(RealVector x) { double xx = x.getEntry(0); double yy = x.getEntry(1); From 7c03f756c08913bfaff2a47c0d0d4aa531757cde Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 10 Jun 2026 08:30:10 -0400 Subject: [PATCH 6/8] use layer and wire alignment in AHDC geometry --- .../detector/alert/AHDC/AlertDCDetector.java | 30 +++ .../detector/alert/AHDC/AlertDCFactory.java | 29 +++ .../alert/AHDC/AlertDCWireIdentifier.java | 244 ++++++++++++++++++ .../main/java/org/jlab/rec/ahdc/Hit/Hit.java | 50 ++-- .../org/jlab/service/ahdc/AHDCEngine.java | 9 +- .../org/jlab/service/alert/ALERTEngine.java | 2 + 6 files changed, 333 insertions(+), 31 deletions(-) create mode 100644 common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCWireIdentifier.java diff --git a/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCDetector.java b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCDetector.java index 14f2220b49..cf570a7f6a 100644 --- a/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCDetector.java +++ b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCDetector.java @@ -7,6 +7,8 @@ import org.jlab.geom.DetectorId; import org.jlab.geom.abs.AbstractDetector; +import org.jlab.geom.prim.Line3D; +import org.jlab.geom.prim.Point3D; /** * @author sergeyeva @@ -26,4 +28,32 @@ public String getType() { return "ALERT DC Detector"; } + public void print() { + System.out.println("/////////////////////////"); + System.out.println("AHDC geometry"); + System.out.println(""); + System.out.println("s : sector"); + System.out.println("sl : super layer"); + System.out.println("l : layer"); + System.out.println("c : component"); + System.out.println("/////////////////////////"); + System.out.println("------------------------------------------------------------------------------"); + System.out.println(" | origin | end"); + System.out.println("------------------------------------------------------------------------------"); + System.out.println("s sl l c | x y z | x y z"); + System.out.println("------------------------------------------------------------------------------"); + for (int s = 1; s <= getNumSectors(); s++) { + for (int sl = 1; sl <= getSector(s).getNumSuperlayers(); sl++) { + for (int l = 1; l <= getSector(s).getSuperlayer(sl).getNumLayers(); l++) { + for (int c = 1; c <= getSector(s).getSuperlayer(sl).getLayer(l).getNumComponents(); c++) { + Line3D line = getSector(s).getSuperlayer(sl).getLayer(l).getComponent(c).getLine(); + Point3D end = line.end(); + Point3D origin = line.origin(); + System.out.printf("%2d %2d %2d %2d | %7.3f %7.3f %7.3f | %7.3f %7.3f %7.3f\n", s, sl, l, c, origin.x(), origin.y(), origin.z(), end.x(), end.y(), end.z()); + } + } + } + } + } + } diff --git a/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java index 3cc32878cd..67c43da27b 100644 --- a/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java +++ b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCFactory.java @@ -205,6 +205,35 @@ public AlertDCLayer createLayer(ConstantProvider cp, int sectorId, int superlaye Point3D p_9 = new Point3D(px_9, py_9, z_end); Point3D p_10 = new Point3D(px_10, py_10, z_end); Point3D p_11 = new Point3D(px_11, py_11, z_end); + + // !!! Alignment correction + int layer_row = AlertDCWireIdentifier.layer2number((superlayerId+1)*10 + (layerId+1)); + double upstream_rotZ = Math.toRadians(cp.getDouble("/geometry/alert/ahdc/layer_alignment/upstream_rotZ", layer_row)); + double downstream_rotZ = Math.toRadians(cp.getDouble("/geometry/alert/ahdc/layer_alignment/downstream_rotZ", layer_row)); + int wire_row = AlertDCWireIdentifier.slc2wire(sectorId+1, (superlayerId+1)*10 + (layerId+1), wireId+1); + double wire_rotZ = Math.toRadians(cp.getDouble("/geometry/alert/ahdc/wire_alignment/rotZ", wire_row)); + + //System.out.printf("%d %d %d, %f , %f, %f\n",sectorId+1, (superlayerId+1)*10 + (layerId+1), wireId+1, Math.toDegrees(upstream_rotZ), Math.toDegrees(downstream_rotZ), Math.toDegrees(wire_rotZ)); + + p_0.rotateZ(wire_rotZ + upstream_rotZ); + p_1.rotateZ(wire_rotZ + upstream_rotZ); + p_2.rotateZ(wire_rotZ + upstream_rotZ); + p_3.rotateZ(wire_rotZ + upstream_rotZ); + p_4.rotateZ(wire_rotZ + upstream_rotZ); + p_5.rotateZ(wire_rotZ + upstream_rotZ); + + p_6.rotateZ(wire_rotZ + downstream_rotZ); + p_7.rotateZ(wire_rotZ + downstream_rotZ); + p_8.rotateZ(wire_rotZ + downstream_rotZ); + p_9.rotateZ(wire_rotZ + downstream_rotZ); + p_10.rotateZ(wire_rotZ + downstream_rotZ); + p_11.rotateZ(wire_rotZ + downstream_rotZ); + + lPoint.rotateZ(upstream_rotZ); + rPoint.rotateZ(downstream_rotZ); + wireLine = new Line3D(lPoint, rPoint); + + // defining a cell around a wireLine, must be counter-clockwise! firstF.add(p_0); firstF.add(p_5); diff --git a/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCWireIdentifier.java b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCWireIdentifier.java new file mode 100644 index 0000000000..fd892eb512 --- /dev/null +++ b/common-tools/clas-geometry/src/main/java/org/jlab/geom/detector/alert/AHDC/AlertDCWireIdentifier.java @@ -0,0 +1,244 @@ +package org.jlab.geom.detector.alert.AHDC; + +public class AlertDCWireIdentifier { + /** Number between 0 and 575 */ + private int num; + + /** Always 1 */ + private int sector; + + /** Can be 11, 21, 22, 31, 32, 41, 42, 51 */ + private int layer; + + /** component id on a given layer, numerotation starting at 1 */ + private int component; + + /** + * Ahdc wire id defined with a num ranging from 0 to 575 + * @param _num + */ + public AlertDCWireIdentifier(int _num) { + num = _num; + int[] res = wire2slc(_num); + sector = res[0]; + layer = res[1]; + component = res[2]; + } + + /** + * Ahdc wire id defined with sector, layer, component identifiers + * @param _sector + * @param _layer can be 11, 21, 22, 31, 32, 41, 51 + * @param _component + */ + public AlertDCWireIdentifier(int _sector, int _layer, int _component) { + sector = _sector; + layer = _layer; + component = _layer; + num = slc2wire(_sector, _layer, _component); + } + + public int getNumber() { return num;} + + public int getSectorId() {return sector;} + + public int getLayerId() {return layer;} + + public int getComponentId() {return component;} + + /** + * Convert wire number (number from 0 to 575) to (sector,layer,component) ids + * + * This is the invert operation of {@link #slc2wire(int, int, int)} + * + * @param wire wire number between 0 and 575 + * @return a triplet (sector, layer, component) in int[] + */ + public static int[] wire2slc(int wire) { + int layer = -1; + int component = -1; + if (wire < 47) { + layer = 11; + component = wire + 1; + } + else if ((47 <= wire) && (wire < 47 + 56)) { + layer = 21; + component = wire - 47 + 1; + } + else if ((47 + 56 <= wire) && (wire < 47 + 56 + 56)) { + layer = 22; + component = wire - 47 - 56 + 1; + } + else if ((47 + 56 + 56 <= wire) && (wire < 47 + 56 + 56 + 72)) { + layer = 31; + component = wire - 47 - 56 - 56 + 1; + } + else if ((47 + 56 + 56 + 72 <= wire) && (wire < 47 + 56 + 56 + 72 + 72)) { + layer = 32; + component = wire - 47 - 56 - 56 - 72 + 1; + } + else if ((47 + 56 + 56 + 72 + 72 <= wire) && (wire < 47 + 56 + 56 + 72 + 72 + 87)) { + layer = 41; + component = wire - 47 - 56 - 56 - 72 - 72 + 1; + } + else if ((47 + 56 + 56 + 72 + 72 + 87 <= wire) && (wire < 47 + 56 + 56 + 72 + 72 + 87 + 87)) { + layer = 42; + component = wire - 47 - 56 - 56 - 72 - 72 - 87 + 1; + } + else { // ((47 + 56 + 56 + 72 + 72 + 87 + 87 <= wire) && (wire < 47 + 56 + 56 + 72 + 72 + 87 + 87 + 99)) { + layer = 51; + component = wire - 47 - 56 - 56 - 72 - 72 - 87 - 87 + 1; + } + return new int[] {1, layer, component}; + } + + /** + * Convert (sector, layer, component) to a unique wire id (number betwwen 0 and 575) + * + * @param sector (not used) + * @param layer + * @param component + * @return unique wire id + */ + public static int slc2wire(int sector, int layer, int component) { + if (layer == 11) { + return component - 1; + } + else if (layer == 21) { + return 47 + component - 1; + } + else if (layer == 22) { + return 47 + 56 + component - 1; + } + else if (layer == 31) { + return 47 + 56 + 56 + component - 1; + } + else if (layer == 32) { + return 47 + 56 + 56 + 72 + component - 1; + } + else if (layer == 41) { + return 47 + 56 + 56 + 72 + 72 + component - 1; + } + else if (layer == 42) { + return 47 + 56 + 56 + 72 + 72 + 87 + component - 1; + } + else if (layer == 51) { + return 47 + 56 + 56 + 72 + 72 + 87 + 87 + component - 1; + } else { + return -1; // not a ahdc wire + } + } + + /** + * Convert the layer digits (11,21,...,51) to layer number between 0 and 7 + * + * @param digit + * @return layer number + */ + public static int layer2number(int digit) { + if (digit == 11) { + return 0; + } + else if (digit == 21) { + return 1; + } + else if (digit == 22) { + return 2; + } + else if (digit == 31) { + return 3; + } + else if (digit == 32) { + return 4; + } + else if (digit == 41) { + return 5; + } + else if (digit == 42) { + return 6; + } + else if (digit == 51) { + return 7; + } else { + return -1; + } + } + + /** + * Convert layer number (from 0 to 7) to the superlayer-layer id (11,21,...,51) + * + * @param digit between 1 and 8 + * @return layer number between (11,21,...,51) + */ + public static int number2layer(int num) { + if (num == 0) { + return 11; + } + else if (num == 1) { + return 21; + } + else if (num == 2) { + return 22; + } + else if (num == 3) { + return 31; + } + else if (num == 4) { + return 32; + } + else if (num == 5) { + return 41; + } + else if (num == 6) { + return 42; + } + else if (num == 7) { + return 51; + } else { + return -1; + } + } + + /** + * + * @param _layer (number 11, 21, 22, ..., 51) + * @return the radius of the _layer + */ + public static double layer2Radius(int _layer) { + if (_layer == 11) { + return 32.0; + } + else if (_layer == 21) { + return 38.0; + } + else if (_layer == 22) { + return 42.0; + } + else if (_layer == 31) { + return 48.0; + } + else if (_layer == 32) { + return 52.0; + } + else if (_layer == 41) { + return 58.0; + } + else if (_layer == 42) { + return 62.0; + } + else if (_layer == 51) { + return 68.0; + } else { + return 0.0; + } + } + + /** + * + * @param _layer_num between 1 and 8 + * @return the radius of the _layer. See {@link #layer2Radius(int)} + */ + public static double layerNum2Radius(int _layer_num) { + return layer2Radius(number2layer(_layer_num)); + } +} 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 cd97e04967..bb1da2f179 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 @@ -250,13 +250,32 @@ public double distance(Point3D point3D) { } public static void main(String[] args) { - AlertDCDetector factory = (new AlertDCFactory()).createDetectorCLAS(new DatabaseConstantProvider()); + int run = 22712; + String variation = "default"; + DatabaseConstantProvider cp = new DatabaseConstantProvider(run, variation); + cp.loadTable("/geometry/alert/ahdc/layer_alignment"); + cp.loadTable("/geometry/alert/ahdc/wire_alignment"); + String upstream_rotZ = "/geometry/alert/ahdc/layer_alignment/upstream_rotZ"; + + System.out.println(cp.getDouble(upstream_rotZ, 0)); + System.out.println(cp.getDouble(upstream_rotZ, 1)); + System.out.println(cp.getDouble(upstream_rotZ, 2)); + System.out.println(cp.getDouble(upstream_rotZ, 3)); + System.out.println(cp.getDouble(upstream_rotZ, 4)); + System.out.println(cp.getDouble(upstream_rotZ, 5)); + System.out.println(cp.getDouble(upstream_rotZ, 6)); + System.out.println(cp.getDouble(upstream_rotZ, 7)); + + + AlertDCDetector AHDCdet = (new AlertDCFactory()).createDetectorCLAS(cp); + //AHDCdet.print(); + 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); + h1.setWirePosition(AHDCdet); + h2.setWirePosition(AHDCdet); System.out.println("h1 : " + h1); System.out.println("h2 : " + h2); System.out.println("h3 : " + h3); @@ -265,31 +284,6 @@ public static void main(String[] args) { System.out.println("h2 compare to h1 : " + h2.compareTo(h1)); System.out.println("h1 compare to h3 : " + h1.compareTo(h3)); - System.out.println("/////////////////////////"); - System.out.println("Test AHDC geometry"); - System.out.println(""); - System.out.println("s : sector"); - System.out.println("sl : super layer"); - System.out.println("l : layer"); - System.out.println("c : component"); - System.out.println("/////////////////////////"); - System.out.println("------------------------------------------------------------------------------"); - System.out.println(" | origin | end"); - System.out.println("------------------------------------------------------------------------------"); - System.out.println("s sl l c | x y z | x y z"); - System.out.println("------------------------------------------------------------------------------"); - for (int s = 1; s <= factory.getNumSectors(); s++) { - for (int sl = 1; sl <= factory.getSector(s).getNumSuperlayers(); sl++) { - for (int l = 1; l <= factory.getSector(s).getSuperlayer(sl).getNumLayers(); l++) { - for (int c = 1; c <= factory.getSector(s).getSuperlayer(sl).getLayer(l).getNumComponents(); c++) { - Line3D line = factory.getSector(s).getSuperlayer(sl).getLayer(l).getComponent(c).getLine(); - Point3D end = line.end(); - Point3D origin = line.origin(); - System.out.printf("%2d %2d %2d %2d | %7.3f %7.3f %7.3f | %7.3f %7.3f %7.3f\n", s, sl, l, c, origin.x(), origin.y(), origin.z(), end.x(), end.y(), end.z()); - } - } - } - } } } 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 992b61919a..5338b146d3 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 @@ -30,7 +30,7 @@ public class AHDCEngine extends ReconstructionEngine { static final Logger LOGGER = Logger.getLogger(AHDCEngine.class.getName()); - private AlertDCDetector factory = null; + private AlertDCDetector AHDCdet = null; private ModeAHDC ahdcExtractor = new ModeAHDC(); // AHDC calibration tables (instance-level, refreshed on run change) @@ -46,7 +46,10 @@ public class AHDCEngine extends ReconstructionEngine { @Override public void detectorChanged(int run) { - factory = (new AlertDCFactory()).createDetectorCLAS(new DatabaseConstantProvider(run,"default")); + DatabaseConstantProvider cp = new DatabaseConstantProvider(run, "default"); + cp.loadTable("/geometry/alert/ahdc/layer_alignment"); + cp.loadTable("/geometry/alert/ahdc/wire_alignment"); + AHDCdet = (new AlertDCFactory()).createDetectorCLAS(cp); } @Override @@ -90,7 +93,7 @@ public boolean processDataEventUser(DataEvent event) { if (event.hasBank("AHDC::adc")) { boolean simulation = event.hasBank("MC::Particle"); - HitReader hitReader = new HitReader(event, factory, simulation, + HitReader hitReader = new HitReader(event, AHDCdet, simulation, ahdcRawHitCutsTable, ahdcTimeOffsetsTable, ahdcTimeToDistanceWireTable, ahdcTimeOverThresholdTable, ahdcAdcGainsTable); ArrayList AHDC_Hits = hitReader.get_AHDCHits(); diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index ea1ee24bde..ad4173cb12 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -141,6 +141,8 @@ public ALERTEngine() { @Override public void detectorChanged(int run) { DatabaseConstantProvider cp = new DatabaseConstantProvider(run, "default"); + cp.loadTable("/geometry/alert/ahdc/layer_alignment"); + cp.loadTable("/geometry/alert/ahdc/wire_alignment"); ATOF = (new AlertTOFFactory()).createDetectorCLAS(cp); AHDC = (new AlertDCFactory()).createDetectorCLAS(cp); } From 70707b3876af196c0c91722678c6dfba272d4690 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 10 Jun 2026 10:04:23 -0400 Subject: [PATCH 7/8] use raw adc in the Kalman Filter --- .../src/main/java/org/jlab/service/alert/ALERTEngine.java | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index ad4173cb12..fc84c9e2e1 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -574,8 +574,9 @@ public boolean processDataEventUser(DataEvent event) { double doca = hitBank.getDouble("doca", hit_row); double time = hitBank.getDouble("time", hit_row); double tot = hitBank.getDouble("timeOverThreshold", hit_row); - // warning : adc is the calibrated one, we need the adc for the Kalman filter - Hit hit = new Hit(id, superlayer, layer, wire, doca, adc, time); + // use raw adc in the Kalman Filter + double gainCorr = ahdcAdcGainsTable.getDoubleValue("gainCorr", 1, 10*superlayer+layer, wire); + Hit hit = new Hit(id, superlayer, layer, wire, doca, adc/gainCorr, time); hit.setWirePosition(AHDC); hit.setTrackId(trackid); hit.setADC(adc); From efb265e5b88a5976e2d2e99c1484dbf57235ceb5 Mon Sep 17 00:00:00 2001 From: Felix Touchte Codjo Date: Wed, 10 Jun 2026 11:14:55 -0400 Subject: [PATCH 8/8] store residual LR in AHDC::hits --- etc/bankdefs/hipo4/alert.json | 4 ++++ .../org/jlab/rec/ahdc/Banks/RecoBankWriter.java | 2 +- .../src/main/java/org/jlab/rec/ahdc/Hit/Hit.java | 9 +++++++++ .../org/jlab/rec/ahdc/KalmanFilter/KFitter.java | 15 +++++++++++++++ .../jlab/rec/ahdc/KalmanFilter/KalmanFilter.java | 1 + 5 files changed, 30 insertions(+), 1 deletion(-) diff --git a/etc/bankdefs/hipo4/alert.json b/etc/bankdefs/hipo4/alert.json index 580a54aecc..96c8b1666c 100644 --- a/etc/bankdefs/hipo4/alert.json +++ b/etc/bankdefs/hipo4/alert.json @@ -228,6 +228,10 @@ "name": "residual", "type": "D", "info": "residual (mm) calculated by the Kalman Filter" + }, { + "name": "residual_LR", + "type": "D", + "info": "residual Left-Right (mm) calculated by the Kalman Filter and used for the AHDC alignment" }, { "name": "time", "type": "D", 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 b276e4b605..82f0461748 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 @@ -12,7 +12,6 @@ import org.jlab.rec.alert.Track.Track; import java.util.ArrayList; -import java.util.List; public class RecoBankWriter { @@ -29,6 +28,7 @@ public DataBank fillAHDCHitsBank(DataEvent event, ArrayList hitList) { bank.setInt("wire", i, hitList.get(i).getWireId()); bank.setDouble("doca", i, hitList.get(i).getDoca()); bank.setDouble("residual", i, hitList.get(i).getResidual()); + bank.setDouble("residual_LR", i, hitList.get(i).getResidual_LR()); bank.setDouble("time", i, hitList.get(i).getTime()); bank.setInt("adc", i, (int) hitList.get(i).getADC()); bank.setDouble("timeOverThreshold", i, hitList.get(i).getToT()); 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 bb1da2f179..3ee813f9b8 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 @@ -32,6 +32,7 @@ public class Hit implements Comparable, KFHit { private double x; private double y; private double residual; + private double residual_LR; private int trackId; //updated constructor with ADC @@ -160,6 +161,14 @@ public void setResidual(double resid) { this.residual = resid; } + public double getResidual_LR() { + return residual_LR; + } + + public void setResidual_LR(double resid) { + this.residual_LR = resid; + } + /** Set calibrated ToT */ public void setToT(double _tot) { this.tot = _tot; 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 12eee314bd..3f64680fdf 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 @@ -7,6 +7,7 @@ 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.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.rec.ahdc.Hit.Hit; @@ -122,6 +123,20 @@ public double residual(Hit hit) { return hit.getDoca()-d; } + // specific to AHDC hits + public double residual_LR(Hit hit) { + Line3D line = hit.getLine().distance( new Point3D( stateEstimation.getEntry(0), stateEstimation.getEntry(1), stateEstimation.getEntry(2) ) ); + //Point3D A = line.origin(); // point on the wire line + Point3D B = line.end(); // point on the track + Point3D C = line.lerpPoint(hit.getDoca()/line.length()); // measurement point on [AB) + + double phiC = Math.atan2(C.y(), C.x()); // between -pi and pi rad + double phiB = Math.atan2(B.y(), B.x()); + //double phiA = Math.atan2(A.y(), A.x()); + + return Math.signum(phiC - phiB)*Math.abs(hit.getDoca()-line.length()); + } + public void ResetErrorCovariance(final RealMatrix initialErrorCovariance){ this.errorCovariance = initialErrorCovariance; } 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 371908ec12..258fe60906 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 @@ -142,6 +142,7 @@ public void propagation(ArrayList tracks, final double magfield, boolean for (Hit hit : AHDC_hits) { PostFitPropagator.predict(hit, true); hit.setResidual(PostFitPropagator.residual(hit)); // output: residual + hit.setResidual_LR(PostFitPropagator.residual_LR(hit)); // output: residual LR } // Fill values for AHDC::kftrack