Skip to content

Commit 9825f60

Browse files
ftouchteftouchte
authored
Move Kalman Filter in ALERT engine (#1100)
* set trackid before the kalman filter * define Niter and the particle type as attributes of the Kalman Filter, make propagation public * remove redundant attributes with _kf extension * move the Kalman Filter in ALERTEngine * add method to set Niter and the particle type * update the AHDC::hits bank, fill residuals * cleaning * uncomment atof part
1 parent 09b42a5 commit 9825f60

File tree

5 files changed

+121
-96
lines changed

5 files changed

+121
-96
lines changed

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Banks/RecoBankWriter.java‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -120,6 +120,7 @@ public DataBank fillAHDCTrackBank(DataEvent event, ArrayList<Track> tracks) {
120120
return bank;
121121
}
122122

123+
/// Here the relavant informations of the tracks are filled in the Kalman Filter
123124
public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> tracks) {
124125

125126
DataBank bank = event.createBank("AHDC::kftrack", tracks.size());
@@ -128,12 +129,12 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> tracks) {
128129

129130
for (Track track : tracks) {
130131
if (track == null) continue;
131-
double x = track.getX0_kf();
132-
double y = track.getY0_kf();
133-
double z = track.getZ0_kf();
134-
double px = track.getPx0_kf();
135-
double py = track.getPy0_kf();
136-
double pz = track.getPz0_kf();
132+
double x = track.get_X0();
133+
double y = track.get_Y0();
134+
double z = track.get_Z0();
135+
double px = track.get_px();
136+
double py = track.get_py();
137+
double pz = track.get_pz();
137138

138139
bank.setInt("trackid", row, (int) track.get_trackId());
139140
bank.setFloat("x", row, (float) x);
@@ -144,9 +145,9 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> tracks) {
144145
bank.setFloat("pz", row, (float) pz);
145146
bank.setInt("n_hits", row, (int) track.get_n_hits());
146147
bank.setInt("sum_adc", row, (int) track.get_sum_adc());
147-
bank.setFloat("path", row, (float) track.get_path_kf());
148-
bank.setFloat("dEdx", row, (float) track.get_dEdx_kf());
149-
bank.setFloat("p_drift", row, (float) track.get_p_drift_kf());
148+
bank.setFloat("path", row, (float) track.get_path());
149+
bank.setFloat("dEdx", row, (float) track.get_dEdx());
150+
bank.setFloat("p_drift", row, (float) track.get_p_drift());
150151
bank.setFloat("chi2", row, (float) track.get_chi2());
151152
bank.setFloat("sum_residuals", row, (float) track.get_sum_residuals());
152153

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/KalmanFilter/KalmanFilter.java‎

Lines changed: 12 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -29,21 +29,23 @@
2929
*/
3030
public class KalmanFilter {
3131

32+
public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;}
33+
3234
public KalmanFilter(ArrayList<Track> tracks, DataEvent event, final double magfield, boolean IsMC) {propagation(tracks, event, magfield, IsMC);}
3335

34-
private final int Niter = 40; // number of iterations for the Kalman Filter
36+
private PDGParticle particle;
37+
private int Niter = 40; // number of iterations for the Kalman Filter
3538
private boolean IsVtxDefined = false; // implemented but not used yet
3639
private double[] vertex_resolutions = {0.09, 1e10}; // {error in r squared in mm^2, error in z squared in mm^2}
3740
// mm, CLAS and AHDC don't necessary have the same alignement (ZERO), this parameter may be subject to calibration
3841
private double clas_alignement = -54;
3942

40-
private void propagation(ArrayList<Track> tracks, DataEvent event, final double magfield, boolean IsMC) {
43+
public void propagation(ArrayList<Track> tracks, DataEvent event, final double magfield, boolean IsMC) {
4144

4245
try {
4346
double vz_constraint = 0; // to be linked to the electron vertex
4447

4548
// Initialization ---------------------------------------------------------------------
46-
final PDGParticle proton = PDGDatabase.getParticleById(2212);
4749
final int numberOfVariables = 6;
4850
final double tesla = 0.001;
4951
final double[] B = {0.0, 0.0, magfield / 10 * tesla};
@@ -74,10 +76,7 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, final double
7476
}
7577

7678
// Loop over tracks
77-
int trackId = 0;
7879
for (Track track : tracks) {
79-
trackId++;
80-
track.set_trackId(trackId);
8180
// Initialize state vector
8281
double x0 = 0.0;
8382
double y0 = 0.0;
@@ -93,7 +92,7 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, final double
9392

9493
// Start propagation
9594
Stepper stepper = new Stepper(y);
96-
RungeKutta4 RK4 = new RungeKutta4(proton, numberOfVariables, B);
95+
RungeKutta4 RK4 = new RungeKutta4(particle, numberOfVariables, B);
9796
Propagator propagator = new Propagator(RK4);
9897

9998
// Initialization of the Kalman Fitter
@@ -139,7 +138,7 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, final double
139138

140139

141140
RealVector x_out = TrackFitter.getStateEstimationVector();
142-
track.setPositionAndMomentumForKF(x_out);
141+
track.setPositionAndMomentumVec(x_out.toArray());
143142

144143
// Post fit propagation (no correction) to set the residuals
145144
KFitter PostFitPropagator = new KFitter(TrackFitter.getStateEstimationVector(), initialErrorCovariance, new Stepper(TrackFitter.getStateEstimationVector().toArray()), new Propagator(RK4), materialHashMap);
@@ -158,23 +157,23 @@ private void propagation(ArrayList<Track> tracks, DataEvent event, final double
158157
double sum_residuals = 0;
159158
double chi2 = 0;
160159
for (Hit hit : AHDC_hits) {
161-
hit.setTrackId(trackId);
162160
sum_adc += hit.getADC();
163161
sum_residuals += hit.getResidual();
164162
chi2 += Math.pow(hit.getResidual(),2)/hit.get_MeasurementNoise().getEntry(0,0);
165163
}
166164
track.set_sum_adc(sum_adc);
167165
track.set_sum_residuals(sum_residuals);
168166
track.set_chi2(chi2/(AHDC_hits.size()-3));
169-
track.set_p_drift_kf(p_drift);
170-
track.set_dEdx_kf(sum_adc/s);
171-
track.set_path_kf(s);
167+
track.set_p_drift(p_drift);
168+
track.set_dEdx(sum_adc/s);
169+
track.set_path(s);
172170
track.set_n_hits(AHDC_hits.size());
173171
}//end of loop on track candidates
174172
} catch (Exception e) {
175173
//e.printStackTrace();
176174
//System.out.println("======> Kalman Filter Error");
177175
}
178176
}
179-
177+
void set_Niter(int Niter) {this.Niter = Niter;}
178+
void set_particle(PDGParticle particle) {this.particle = particle;}
180179
}

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java‎

Lines changed: 5 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -35,16 +35,7 @@ public class Track {
3535
private double dEdx = 0; ///< deposited energy per path length (adc/mm)
3636
private double p_drift = 0; ///< momentum in the drift region (MeV)
3737
private double path = 0; ///< length of the track (mm)
38-
// AHDC::kftrack
39-
private double x0_kf = 0;
40-
private double y0_kf = 0;
41-
private double z0_kf = 0;
42-
private double px0_kf = 0;
43-
private double py0_kf = 0;
44-
private double pz0_kf = 0;
45-
private double dEdx_kf = 0; ///< deposited energy per path length (adc/mm)
46-
private double p_drift_kf = 0; ///< momentum in the drift region (MeV)
47-
private double path_kf = 0; ///< length of the track (mm)
38+
4839
// AHDC::aiprediction
4940
private int predicted_ATOF_sector = -1;
5041
private int predicted_ATOF_layer = -1;
@@ -93,15 +84,6 @@ public void setPositionAndMomentumVec(double[] x) {
9384
this.pz0 = x[5];
9485
}
9586

96-
public void setPositionAndMomentumForKF(RealVector x) {
97-
this.x0_kf = x.getEntry(0);
98-
this.y0_kf = x.getEntry(1);
99-
this.z0_kf = x.getEntry(2);
100-
this.px0_kf = x.getEntry(3);
101-
this.py0_kf = x.getEntry(4);
102-
this.pz0_kf = x.getEntry(5);
103-
}
104-
10587
private void generateHitList() {
10688
for (Cluster cluster : _Clusters) {
10789
for (PreCluster preCluster : cluster.get_PreClusters_list()) {
@@ -174,31 +156,6 @@ public double get_pz() {
174156
return pz0;
175157
}
176158

177-
public double getX0_kf() {
178-
return x0_kf;
179-
}
180-
181-
public double getY0_kf() {
182-
return y0_kf;
183-
}
184-
185-
public double getZ0_kf() {
186-
return z0_kf;
187-
}
188-
189-
public double getPx0_kf() {
190-
return px0_kf;
191-
}
192-
193-
public double getPy0_kf() {
194-
return py0_kf;
195-
}
196-
197-
public double getPz0_kf() {
198-
return pz0_kf;
199-
}
200-
201-
// Same for Track and KFTrack
202159
public void set_trackId(int _trackId) {
203160
trackId = _trackId;
204161
// set trackId for clusters
@@ -209,6 +166,10 @@ public void set_trackId(int _trackId) {
209166
for(InterCluster interCluster : this._InterClusters) {
210167
interCluster.setTrackId(_trackId);
211168
}
169+
// set trackId for hits
170+
for (Hit hit : this.hits) {
171+
hit.setTrackId(_trackId);
172+
}
212173
}
213174
public void set_n_hits(int _n_hits) { n_hits = _n_hits;}
214175
public void set_sum_adc(int _sum_adc) { sum_adc = _sum_adc;}
@@ -226,14 +187,6 @@ public void set_trackId(int _trackId) {
226187
public double get_dEdx() {return dEdx;}
227188
public double get_p_drift() {return p_drift;}
228189
public double get_path() {return path;}
229-
230-
// AHDC::kftrack
231-
public void set_dEdx_kf(double _dEdx_kf) { dEdx_kf = _dEdx_kf;}
232-
public void set_p_drift_kf(double _p_drift_kf) { p_drift_kf = _p_drift_kf;}
233-
public void set_path_kf(double _path_kf) { path_kf = _path_kf;}
234-
public double get_dEdx_kf() {return dEdx_kf;}
235-
public double get_p_drift_kf() {return p_drift_kf;}
236-
public double get_path_kf() {return path_kf;}
237190

238191
// AHDC::aiprediction
239192
public void set_predicted_ATOF_sector(int s) {predicted_ATOF_sector = s;}

‎reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java‎

Lines changed: 5 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,6 @@
1515
import org.jlab.rec.ahdc.Hit.Hit;
1616
import org.jlab.rec.ahdc.Hit.HitReader;
1717
import org.jlab.rec.ahdc.HoughTransform.HoughTransform;
18-
import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter;
1918
import org.jlab.rec.ahdc.KalmanFilter.MaterialMap;
2019
import org.jlab.rec.ahdc.PreCluster.PreCluster;
2120
import org.jlab.rec.ahdc.PreCluster.PreClusterFinder;
@@ -100,7 +99,7 @@ else if (Objects.equals(this.getEngineConfigString("Mode"), ModeTrackFinding.CV_
10099

101100
this.getConstantsManager().setVariation("default");
102101

103-
this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction");
102+
this.registerOutputBank("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::mc","AHDC::ai:prediction");
104103

105104
return true;
106105
}
@@ -110,16 +109,13 @@ else if (Objects.equals(this.getEngineConfigString("Mode"), ModeTrackFinding.CV_
110109
@Override
111110
public boolean processDataEvent(DataEvent event) {
112111

113-
double magfield = 50.0; // what is this? The full magnetic field strength in kGauss (factor * 50kGauss)
114-
115112
if(event.hasBank("MC::Particle")) simulation = true;
116113

117114
ahdcExtractor.update(30, null, event, "AHDC::wf", "AHDC::adc");
118115

119116
if (event.hasBank("RUN::config")) {
120117
DataBank bank = event.getBank("RUN::config");
121118
int newRun = bank.getInt("run", 0);
122-
float magfieldfactor = bank.getFloat("solenoid", 0);
123119
if (newRun <= 0) {
124120
LOGGER.warning("AHDCEngine: got run <= 0 in RUN::config, skipping event.");
125121
return false;
@@ -130,10 +126,6 @@ public boolean processDataEvent(DataEvent event) {
130126
CalibrationConstantsLoader.Load(newRun, this.getConstantsManager());
131127
Run = newRun;
132128
}
133-
134-
/// What is this? The field value in the RUN::config bank is a scaling factor (between -1 and 1) of the full field
135-
/// The kalman filter use the field in kG not Tesla
136-
magfield = 50 * magfieldfactor;
137129
}
138130

139131

@@ -222,7 +214,10 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) {
222214
//AHDC_Tracks.add(new Track(AHDC_Hits));
223215

224216
// V) Global fit
217+
int trackid = 0;
225218
for (Track track : AHDC_Tracks) {
219+
trackid++;
220+
track.set_trackId(trackid);
226221
int nbOfPoints = track.get_Clusters().size();
227222

228223
double[][] szPos = new double[nbOfPoints][3];
@@ -239,9 +234,6 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) {
239234
track.setPositionAndMomentum(h.HelixFit(nbOfPoints, szPos, 1));
240235
}
241236

242-
// VI) Kalman Filter
243-
KalmanFilter kalmanFitter = new KalmanFilter(AHDC_Tracks, event, magfield, simulation);
244-
245237
// VII) Write bank
246238
RecoBankWriter writer = new RecoBankWriter();
247239

@@ -253,7 +245,6 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) {
253245
}
254246
DataBank recoClusterBank = writer.fillClustersBank(event, AHDC_Clusters);
255247
DataBank recoTracksBank = writer.fillAHDCTrackBank(event, AHDC_Tracks);
256-
DataBank recoKFTracksBank = writer.fillAHDCKFTrackBank(event, AHDC_Tracks);
257248

258249
ArrayList<InterCluster> all_interclusters = new ArrayList<>();
259250
for (Track track : AHDC_Tracks) {
@@ -262,11 +253,11 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) {
262253
DataBank recoInterClusterBank = writer.fillInterClusterBank(event, all_interclusters);
263254
// DataBank AIPredictionBanks = writer.fillAIPrediction(event, predictions);
264255

256+
//event.removeBanks("AHDC::hits","AHDC::preclusters","AHDC::clusters","AHDC::track","AHDC::kftrack","AHDC::mc","AHDC::ai:prediction");
265257
event.appendBank(recoHitsBank);
266258
event.appendBank(recoPreClusterBank);
267259
event.appendBank(recoClusterBank);
268260
event.appendBank(recoTracksBank);
269-
event.appendBank(recoKFTracksBank);
270261
event.appendBank(recoInterClusterBank);
271262
// event.appendBank(AIPredictionBanks);
272263

0 commit comments

Comments
 (0)