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 6adf5a1ce4..5c32ecf73b 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 @@ -120,6 +120,7 @@ public DataBank fillAHDCTrackBank(DataEvent event, ArrayList tracks) { return bank; } + /// Here the relavant informations of the tracks are filled in the Kalman Filter public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { DataBank bank = event.createBank("AHDC::kftrack", tracks.size()); @@ -128,12 +129,12 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { for (Track track : tracks) { if (track == null) continue; - double x = track.getX0_kf(); - double y = track.getY0_kf(); - double z = track.getZ0_kf(); - double px = track.getPx0_kf(); - double py = track.getPy0_kf(); - double pz = track.getPz0_kf(); + double x = track.get_X0(); + double y = track.get_Y0(); + double z = track.get_Z0(); + double px = track.get_px(); + double py = track.get_py(); + double pz = track.get_pz(); bank.setInt("trackid", row, (int) track.get_trackId()); bank.setFloat("x", row, (float) x); @@ -144,9 +145,9 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList tracks) { bank.setFloat("pz", row, (float) pz); bank.setInt("n_hits", row, (int) track.get_n_hits()); bank.setInt("sum_adc", row, (int) track.get_sum_adc()); - bank.setFloat("path", row, (float) track.get_path_kf()); - bank.setFloat("dEdx", row, (float) track.get_dEdx_kf()); - bank.setFloat("p_drift", row, (float) track.get_p_drift_kf()); + bank.setFloat("path", row, (float) track.get_path()); + bank.setFloat("dEdx", row, (float) track.get_dEdx()); + bank.setFloat("p_drift", row, (float) track.get_p_drift()); bank.setFloat("chi2", row, (float) track.get_chi2()); bank.setFloat("sum_residuals", row, (float) track.get_sum_residuals()); 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 4d2840ee73..272857efbd 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 @@ -29,21 +29,23 @@ */ public class KalmanFilter { + public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;} + public KalmanFilter(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);} - private final int Niter = 40; // number of iterations for the Kalman Filter + 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 = -54; - private void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { + public void propagation(ArrayList tracks, DataEvent event, final double magfield, boolean IsMC) { try { double vz_constraint = 0; // to be linked to the electron vertex // Initialization --------------------------------------------------------------------- - final PDGParticle proton = PDGDatabase.getParticleById(2212); final int numberOfVariables = 6; final double tesla = 0.001; final double[] B = {0.0, 0.0, magfield / 10 * tesla}; @@ -74,10 +76,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double } // Loop over tracks - int trackId = 0; for (Track track : tracks) { - trackId++; - track.set_trackId(trackId); // Initialize state vector double x0 = 0.0; double y0 = 0.0; @@ -93,7 +92,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double // Start propagation Stepper stepper = new Stepper(y); - RungeKutta4 RK4 = new RungeKutta4(proton, numberOfVariables, B); + RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B); Propagator propagator = new Propagator(RK4); // Initialization of the Kalman Fitter @@ -139,7 +138,7 @@ private void propagation(ArrayList tracks, DataEvent event, final double RealVector x_out = TrackFitter.getStateEstimationVector(); - track.setPositionAndMomentumForKF(x_out); + track.setPositionAndMomentumVec(x_out.toArray()); // 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); @@ -158,7 +157,6 @@ private void propagation(ArrayList tracks, DataEvent event, final double double sum_residuals = 0; double chi2 = 0; for (Hit hit : AHDC_hits) { - hit.setTrackId(trackId); sum_adc += hit.getADC(); sum_residuals += hit.getResidual(); chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0); @@ -166,9 +164,9 @@ private void propagation(ArrayList tracks, DataEvent event, final double track.set_sum_adc(sum_adc); track.set_sum_residuals(sum_residuals); track.set_chi2(chi2/(AHDC_hits.size()-3)); - track.set_p_drift_kf(p_drift); - track.set_dEdx_kf(sum_adc/s); - track.set_path_kf(s); + track.set_p_drift(p_drift); + track.set_dEdx(sum_adc/s); + track.set_path(s); track.set_n_hits(AHDC_hits.size()); }//end of loop on track candidates } catch (Exception e) { @@ -176,5 +174,6 @@ private void propagation(ArrayList tracks, DataEvent event, final double //System.out.println("======> Kalman Filter Error"); } } - + void set_Niter(int Niter) {this.Niter = Niter;} + void set_particle(PDGParticle particle) {this.particle = particle;} } diff --git a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java index 5d140c3c45..e9944b7846 100644 --- a/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java +++ b/reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java @@ -35,16 +35,7 @@ public class Track { private double dEdx = 0; ///< deposited energy per path length (adc/mm) private double p_drift = 0; ///< momentum in the drift region (MeV) private double path = 0; ///< length of the track (mm) - // AHDC::kftrack - private double x0_kf = 0; - private double y0_kf = 0; - private double z0_kf = 0; - private double px0_kf = 0; - private double py0_kf = 0; - private double pz0_kf = 0; - private double dEdx_kf = 0; ///< deposited energy per path length (adc/mm) - private double p_drift_kf = 0; ///< momentum in the drift region (MeV) - private double path_kf = 0; ///< length of the track (mm) + // AHDC::aiprediction private int predicted_ATOF_sector = -1; private int predicted_ATOF_layer = -1; @@ -93,15 +84,6 @@ public void setPositionAndMomentumVec(double[] x) { this.pz0 = x[5]; } - public void setPositionAndMomentumForKF(RealVector x) { - this.x0_kf = x.getEntry(0); - this.y0_kf = x.getEntry(1); - this.z0_kf = x.getEntry(2); - this.px0_kf = x.getEntry(3); - this.py0_kf = x.getEntry(4); - this.pz0_kf = x.getEntry(5); - } - private void generateHitList() { for (Cluster cluster : _Clusters) { for (PreCluster preCluster : cluster.get_PreClusters_list()) { @@ -174,31 +156,6 @@ public double get_pz() { return pz0; } - public double getX0_kf() { - return x0_kf; - } - - public double getY0_kf() { - return y0_kf; - } - - public double getZ0_kf() { - return z0_kf; - } - - public double getPx0_kf() { - return px0_kf; - } - - public double getPy0_kf() { - return py0_kf; - } - - public double getPz0_kf() { - return pz0_kf; - } - - // Same for Track and KFTrack public void set_trackId(int _trackId) { trackId = _trackId; // set trackId for clusters @@ -209,6 +166,10 @@ public void set_trackId(int _trackId) { for(InterCluster interCluster : this._InterClusters) { interCluster.setTrackId(_trackId); } + // set trackId for hits + for (Hit hit : this.hits) { + hit.setTrackId(_trackId); + } } public void set_n_hits(int _n_hits) { n_hits = _n_hits;} public void set_sum_adc(int _sum_adc) { sum_adc = _sum_adc;} @@ -226,14 +187,6 @@ public void set_trackId(int _trackId) { public double get_dEdx() {return dEdx;} public double get_p_drift() {return p_drift;} public double get_path() {return path;} - - // AHDC::kftrack - public void set_dEdx_kf(double _dEdx_kf) { dEdx_kf = _dEdx_kf;} - public void set_p_drift_kf(double _p_drift_kf) { p_drift_kf = _p_drift_kf;} - public void set_path_kf(double _path_kf) { path_kf = _path_kf;} - public double get_dEdx_kf() {return dEdx_kf;} - public double get_p_drift_kf() {return p_drift_kf;} - public double get_path_kf() {return path_kf;} // AHDC::aiprediction public void set_predicted_ATOF_sector(int s) {predicted_ATOF_sector = s;} 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 b7cbb1b01c..6e55b25c88 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 @@ -15,7 +15,6 @@ import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Hit.HitReader; import org.jlab.rec.ahdc.HoughTransform.HoughTransform; -import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; import org.jlab.rec.ahdc.KalmanFilter.MaterialMap; import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.ahdc.PreCluster.PreClusterFinder; @@ -100,7 +99,7 @@ else if (Objects.equals(this.getEngineConfigString("Mode"), ModeTrackFinding.CV_ this.getConstantsManager().setVariation("default"); - this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction"); + this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::mc","AHDC::ai:prediction"); return true; } @@ -110,8 +109,6 @@ else if (Objects.equals(this.getEngineConfigString("Mode"), ModeTrackFinding.CV_ @Override public boolean processDataEvent(DataEvent event) { - double magfield = 50.0; // what is this? The full magnetic field strength in kGauss (factor * 50kGauss) - if(event.hasBank("MC::Particle")) simulation = true; ahdcExtractor.update(30, null, event, "AHDC::wf", "AHDC::adc"); @@ -119,7 +116,6 @@ public boolean processDataEvent(DataEvent event) { if (event.hasBank("RUN::config")) { DataBank bank = event.getBank("RUN::config"); int newRun = bank.getInt("run", 0); - float magfieldfactor = bank.getFloat("solenoid", 0); if (newRun <= 0) { LOGGER.warning("AHDCEngine: got run <= 0 in RUN::config, skipping event."); return false; @@ -130,10 +126,6 @@ public boolean processDataEvent(DataEvent event) { CalibrationConstantsLoader.Load(newRun, this.getConstantsManager()); Run = newRun; } - - /// What is this? The field value in the RUN::config bank is a scaling factor (between -1 and 1) of the full field - /// The kalman filter use the field in kG not Tesla - magfield = 50 * magfieldfactor; } @@ -222,7 +214,10 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) { //AHDC_Tracks.add(new Track(AHDC_Hits)); // V) Global fit + int trackid = 0; for (Track track : AHDC_Tracks) { + trackid++; + track.set_trackId(trackid); int nbOfPoints = track.get_Clusters().size(); double[][] szPos = new double[nbOfPoints][3]; @@ -239,9 +234,6 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) { track.setPositionAndMomentum(h.HelixFit(nbOfPoints, szPos, 1)); } - // VI) Kalman Filter - KalmanFilter kalmanFitter = new KalmanFilter(AHDC_Tracks, event, magfield, simulation); - // VII) Write bank RecoBankWriter writer = new RecoBankWriter(); @@ -253,7 +245,6 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) { } DataBank recoClusterBank = writer.fillClustersBank(event, AHDC_Clusters); DataBank recoTracksBank = writer.fillAHDCTrackBank(event, AHDC_Tracks); - DataBank recoKFTracksBank = writer.fillAHDCKFTrackBank(event, AHDC_Tracks); ArrayList all_interclusters = new ArrayList<>(); for (Track track : AHDC_Tracks) { @@ -262,11 +253,11 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) { DataBank recoInterClusterBank = writer.fillInterClusterBank(event, all_interclusters); // DataBank AIPredictionBanks = writer.fillAIPrediction(event, predictions); + //event.removeBanks("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction"); event.appendBank(recoHitsBank); event.appendBank(recoPreClusterBank); event.appendBank(recoClusterBank); event.appendBank(recoTracksBank); - event.appendBank(recoKFTracksBank); event.appendBank(recoInterClusterBank); // event.appendBank(AIPredictionBanks); 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 a5cdc2518e..caeeba638f 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 @@ -1,9 +1,11 @@ package org.jlab.service.alert; +import java.io.File; import java.util.ArrayList; import java.util.concurrent.atomic.AtomicInteger; -import java.io.File; +import org.jlab.clas.reco.ReconstructionEngine; +import org.jlab.clas.swimtools.Swim; import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.base.Detector; import org.jlab.geom.detector.alert.ATOF.AlertTOFFactory; @@ -11,15 +13,19 @@ import org.jlab.io.base.DataEvent; import org.jlab.io.hipo.HipoDataSource; import org.jlab.io.hipo.HipoDataSync; - -import org.jlab.clas.reco.ReconstructionEngine; -import org.jlab.clas.swimtools.Swim; - import org.jlab.rec.alert.TrackMatchingAI.ModelTrackMatching; - import org.jlab.rec.alert.banks.RecoBankWriter; import org.jlab.rec.alert.projections.TrackProjector; import org.jlab.rec.atof.hit.ATOFHit; +import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; +import org.jlab.rec.ahdc.Hit.Hit; +import org.jlab.geom.detector.alert.AHDC.AlertDCDetector; +import org.jlab.geom.detector.alert.AHDC.AlertDCFactory; +import org.jlab.rec.ahdc.Track.Track; +import org.jlab.clas.pdg.PDGDatabase; +import org.jlab.clas.pdg.PDGParticle; + + import ai.djl.util.Pair; @@ -47,6 +53,7 @@ public class ALERTEngine extends ReconstructionEngine { private RecoBankWriter rbc; Detector ATOF; // ALERT ATOF detector + private AlertDCDetector AHDC; // ALERT AHDC detector /** * Current run number being processed. @@ -87,6 +94,7 @@ public boolean init() { AlertTOFFactory factory = new AlertTOFFactory(); DatabaseConstantProvider cp = new DatabaseConstantProvider(11, "default"); ATOF = factory.createDetectorCLAS(cp); + AHDC = (new AlertDCFactory()).createDetectorCLAS(new DatabaseConstantProvider()); if(this.getEngineConfigString("Mode")!=null) { //if (Objects.equals(this.getEngineConfigString("Mode"), Mode.AI_Track_Finding.name())) @@ -116,9 +124,9 @@ public boolean processDataEvent(DataEvent event) { return true; } - DataBank bank = event.getBank("RUN::config"); + DataBank runBank = event.getBank("RUN::config"); - int newRun = bank.getInt("run", 0); + int newRun = runBank.getInt("run", 0); if (newRun == 0) { return true; } @@ -214,6 +222,79 @@ public boolean processDataEvent(DataEvent event) { } } rbc.appendTrackMatchingAIBank(event, matched_ATOF_hit_id); + + /////////////////////////////////////////// + /// Kalmam Filter + /// /////////////////////////////////////// + + // read the list of tracks/hits from the banks AHDC::track and AHDC::hits + if (!event.hasBank("AHDC::track")) {return false;} + DataBank trackBank = event.getBank("AHDC::track"); + DataBank hitBank = event.getBank("AHDC::hits"); + ArrayList AHDC_tracks = new ArrayList<>(); + for (int row = 0; row < trackBank.rows(); row++) { + int trackid = trackBank.getInt("trackid", row); + ArrayList AHDC_hits = new ArrayList<>(); + for (int hit_row = 0; hit_row < hitBank.rows(); hit_row++) { + if(trackid == hitBank.getInt("trackid", hit_row)) { + int id = hitBank.getShort("id", hit_row); + int superlayer = hitBank.getByte("superlayer", hit_row); + int layer = hitBank.getByte("layer", hit_row); + int wire = hitBank.getInt("wire", hit_row); + int adc = hitBank.getInt("adc", hit_row); + 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); + hit.setWirePosition(AHDC); + hit.setTrackId(trackid); + hit.setADC(adc); + hit.setToT(tot); + AHDC_hits.add(hit); + } + } + AHDC_tracks.add(new Track(AHDC_hits)); + // Initialise the position and the momentum using the information of the AHDC::track + // position : mm + // momentum : MeV + double x = trackBank.getFloat("x", row); + double y = trackBank.getFloat("y", row); + double z = trackBank.getFloat("z", row); + double px = trackBank.getFloat("px", row); + double py = trackBank.getFloat("py", row); + double pz = trackBank.getFloat("pz", row); + double[] vec = {x, y, z, px, py, pz}; + AHDC_tracks.get(row).setPositionAndMomentumVec(vec); + AHDC_tracks.get(row).set_trackId(trackid); + } + // 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; + 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 + ///////////////////////////////////////////// + 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 + event.removeBank("AHDC::hits"); + ArrayList AHDC_hits = new ArrayList<>(); + for (Track track : AHDC_tracks) { + AHDC_hits.addAll(track.getHits()); + } + DataBank recoKFHitsBank = ahdc_writer.fillAHDCHitsBank(event, AHDC_hits); + event.appendBank(recoKFHitsBank); // remark: only hits assocuated to a track are saved + + return true; }