Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ public DataBank fillAHDCTrackBank(DataEvent event, ArrayList<Track> tracks) {
return bank;
}

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

DataBank bank = event.createBank("AHDC::kftrack", tracks.size());
Expand All @@ -128,12 +129,12 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> 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);
Expand All @@ -144,9 +145,9 @@ public DataBank fillAHDCKFTrackBank(DataEvent event, ArrayList<Track> 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());

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -29,21 +29,23 @@
*/
public class KalmanFilter {

public KalmanFilter(PDGParticle particle, int Niter) {this.particle = particle; this.Niter = Niter;}

public KalmanFilter(ArrayList<Track> 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<Track> tracks, DataEvent event, final double magfield, boolean IsMC) {
public void propagation(ArrayList<Track> 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};
Expand Down Expand Up @@ -74,10 +76,7 @@ private void propagation(ArrayList<Track> 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;
Expand All @@ -93,7 +92,7 @@ private void propagation(ArrayList<Track> 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
Expand Down Expand Up @@ -139,7 +138,7 @@ private void propagation(ArrayList<Track> 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);
Expand All @@ -158,23 +157,23 @@ private void propagation(ArrayList<Track> 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);
}
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) {
//e.printStackTrace();
//System.out.println("======> Kalman Filter Error");
}
}

void set_Niter(int Niter) {this.Niter = Niter;}
void set_particle(PDGParticle particle) {this.particle = particle;}
}
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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()) {
Expand Down Expand Up @@ -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
Expand All @@ -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;}
Expand All @@ -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;}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand All @@ -110,16 +109,13 @@ 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");

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;
Expand All @@ -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;
}


Expand Down Expand Up @@ -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];
Expand All @@ -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();

Expand All @@ -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<InterCluster> all_interclusters = new ArrayList<>();
for (Track track : AHDC_Tracks) {
Expand All @@ -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);

Expand Down
Loading