Skip to content

Kalman dev #850

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pbutti merged 4 commits into from
May 1, 2021
Merged

Kalman dev #850

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d8b0eca
Modifications to debug plots.
robertprestonjohnson Apr 14, 2021
49f1c75
Made the edge tolerance dependent on the track extrapolation error. F…
robertprestonjohnson Apr 15, 2021
32571c0
Just added some debug plots to study events with bad covariance matri…
robertprestonjohnson Apr 22, 2021
8fd841a
Updated the defaults in KalmanParams to match the steering file used …
robertprestonjohnson Apr 27, 2021
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
4 changes: 2 additions & 2 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,8 +93,8 @@ String toString(String s) {
String str;
str = String.format("HelixState %s: helix parameters=%s, pivot=%s\n", s, a.toString(), X0.toString());
str = str + String.format(" Origin=%s, B=%10.6f in direction %s\n", origin.toString(), B, tB.toString());
str = str + "Covariance:" + C.toString();
str = str + Rot.toString("from global coordinates to field coordinates");
if (C != null) str = str + "Covariance:" + C.toString();
if (Rot != null) str = str + Rot.toString("from global coordinates to field coordinates");
str = str + "End of HelixState dump\n";
return str;
}
Expand Down
15 changes: 10 additions & 5 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code
// Control parameters
// Units are Tesla, GeV, mm

int nTrials = 10000; // The number of test events to generate for fitting
int nTrials = 100; // The number of test events to generate for fitting
int startLayer = 10; // Where to start the Kalman filtering
int nIteration = 2; // Number of filter iterations
int nAxial = 3; // Number of axial layers needed by the linear fit
Expand All @@ -49,6 +49,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code

boolean verbose = nTrials < 2;
double executionTime = 0.;
int nBadCov = 0;

double eCalLoc = 1394.;

Expand Down Expand Up @@ -97,9 +98,9 @@ class HelixTest3 { // Program for testing the Kalman fitting code
Vec B = new Vec(3, fM.getField(new Vec(0., y, z)));
System.out.format("x=0 y=%6.1f z=%6.1f: %s\n", y, z, B.toString());
}
System.out.format("B field map vs x at ECAL:\n");
for (double x=-200.; x<200.; x+=5.) {
double y=eCalLoc + 10.;
System.out.format("B field map vs x at center:\n");
for (double x=-300.; x<300.; x+=5.) {
double y=505.;
double z=20.;
Vec B = new Vec(3, fM.getField(new Vec(x, y, z)));
System.out.format("x=%6.1f y=%6.1f z=%6.1f: %s\n", x, y, z, B.toString());
Expand Down Expand Up @@ -418,7 +419,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code
hResidS4[i] = new Histogram(100, -0.1, 0.002, String.format("Smoothed true residual for plane %d", i), "mm", "hits");
hResidX[i] = new Histogram(100, -0.8, 0.016, String.format("True residual in global X for plane %d", i), "mm", "hits");
hResidZ[i] = new Histogram(100, -0.1, 0.002, String.format("True residual in global Z for plane %d", i), "mm", "hits");
hUnbias[i] = new Histogram(100, -0.2, 0.004, String.format("Unbiased residual for plant %d", i), "mm", "hits");
hUnbias[i] = new Histogram(100, -0.2, 0.004, String.format("Unbiased residual for plane %d", i), "mm", "hits");
hUnbiasSig[i] = new Histogram(100, -10., 0.2, String.format("Unbiased residuals for layer %d", i), "sigmas", "hits");
}
Histogram hpropx = new Histogram(100,-5.,0.1,"projected track-state x error","mm","track");
Expand Down Expand Up @@ -740,11 +741,14 @@ class HelixTest3 { // Program for testing the Kalman fitting code
Matrix C = new Matrix(KalmanTrack.originCovariance());
EigenvalueDecomposition eED= new EigenvalueDecomposition(C);
double [] ev = eED.getRealEigenvalues();
boolean badCov = false;
for (int i=0; i<5; ++i) {
if (ev[i] < 0.) {
System.out.format("Event %d, eigenvalue %d of covariance is negative!", iTrial, i);
badCov = true;
}
}
if (badCov) nBadCov++;
if (iTrial < 10) {
Vec evV = new Vec(5,ev);
evV.print("Eigenvalues of covariance");
Expand Down Expand Up @@ -1107,6 +1111,7 @@ else if (kF.sites.indexOf(site) == kF.sites.size()-1) {
}
}

System.out.format("Number of track fits with bad covariance matrix = %d\n", nBadCov);
long grdef = rnd.nextLong();
System.out.format("New seed = %d\n", grdef);
timestamp = Instant.now();
Expand Down
4 changes: 2 additions & 2 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,8 +100,8 @@ public class KalTrack {
for (int idx=firstSite; idx<=lastSite; ++idx) {
MeasurementSite site = SiteList.get(idx);
if (site.aS == null) { // This should never happen
logger.log(Level.WARNING, String.format("Event %d: site is missing smoothed state vector for layer %d detector %d",
eventNumber, site.m.Layer, site.m.detector));
logger.log(Level.SEVERE, String.format("Event %d: site of track %d is missing smoothed state vector for layer %d detector %d",
eventNumber, ID, site.m.Layer, site.m.detector));
logger.log(Level.WARNING, site.toString("bad site"));
bad = true;
continue;
Expand Down
4 changes: 3 additions & 1 deletion tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ static Vec getField(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) {

static double [] getFielD(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) {
// Field map for stand-alone running
if (FieldMap.class.isInstance(hpsFm)) { return ((FieldMap) (hpsFm)).getField(kalPos); }
if (FieldMap.class.isInstance(hpsFm)) return ((FieldMap) (hpsFm)).getField(kalPos);

// Standard field map for running in hps-java
//System.out.format("Accessing HPS field map for position %8.3f %8.3f %8.3f\n", kalPos.v[0], kalPos.v[1], kalPos.v[2]);
Expand All @@ -130,6 +130,8 @@ static Vec getField(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) {
} else {
if (hpsPos[1] > 70.0) hpsPos[1] = 70.0; // To avoid getting a field returned that is identically equal to zero
if (hpsPos[1] < -70.0) hpsPos[1] = -70.0;
if (hpsPos[0] < -225.) hpsPos[0] = -225.;
if (hpsPos[0] > 270.) hpsPos[0] = 270.;
}
double[] hpsField = hpsFm.getField(hpsPos);
if (uniformB) {
Expand Down
76 changes: 23 additions & 53 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanParams.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ public void print() {
mxTdif = 30.; // Maximum time difference of hits in a track
firstLayer = 0; // First layer in the tracking system (2 for pre-2019 data)
lowPhThresh = 0.25; // Residual improvement ratio necessary to use a low-ph hit instead of high-ph
seedCompThr = -1.; // Remove SeedTracks with all Helix params within relative seedCompThr . If -1 do not apply duplicate removal
seedCompThr = 0.05; // Remove SeedTracks with all Helix params within relative seedCompThr . If -1 do not apply duplicate removal

// Load the default search strategies
// Index 0 is for the bottom tracker (+z), 1 for the top (-z)
Expand All @@ -169,42 +169,24 @@ public void print() {
// A S A S A S A S A S A S A S top
// 0,1,2,3,4,5,6,7,8,9,10,11,12,13
// S A S A S A S A S A S A S A bottom
final int[] list0 = {6, 7, 8, 9, 10};
final int[] list1 = {4, 5, 6, 7, 8};
final int[] list2 = {5, 6, 8, 9, 10};
final int[] list3 = {5, 6, 7, 8, 10};
final int[] list4 = { 3, 6, 8, 9, 10 };
final int[] list5 = { 4, 5, 8, 9, 10 };
final int[] list6 = { 4, 6, 7, 8, 9 };
final int[] list7 = { 4, 6, 7, 9, 10 };
final int[] list8 = { 2, 5, 8, 9, 12};
final int[] list9 = { 8, 10, 11, 12, 13};
final int[] list10 = {6, 9, 10, 11, 12};
final int[] list11 = {6, 7, 9, 10, 12};
final int[] list12 = {2, 3, 4, 5, 6};
final int[] list13 = {2, 4, 5, 6, 7};
final int[] list14 = {6, 7, 8, 10, 11};
final int[] list15 = {1, 2, 3, 4, 6};
final int[] list16 = {0, 2, 3, 4, 5};
final int[] list17 = {0, 3, 4, 5, 6};
lyrList[0].add(list0);
lyrList[0].add(list1);
lyrList[0].add(list2);
lyrList[0].add(list3);
lyrList[0].add(list4);
lyrList[0].add(list5);
lyrList[0].add(list6);
lyrList[0].add(list7);
lyrList[0].add(list8);
lyrList[0].add(list9);
lyrList[0].add(list10);
lyrList[0].add(list11);
lyrList[0].add(list12);
lyrList[0].add(list13);
lyrList[0].add(list14);
lyrList[0].add(list15);
lyrList[0].add(list16);
lyrList[0].add(list17);
addStrategy("000BBS0");
addStrategy("00BBS00");
addStrategy("00ASBS0");
addStrategy("00ABSS0");
addStrategy("0A0SBS0");
addStrategy("00B0BS0");
addStrategy("00SBSA0");
addStrategy("00SBB00");
addStrategy("00SBAS0");
addStrategy("0SA0BS0");
addStrategy("0000SBB");
addStrategy("000SABS");
addStrategy("0BBS000");
addStrategy("0SBB000");
addStrategy("000BSSA");
addStrategy("ABSS000");
addStrategy("SBB0000");
addStrategy("SABS000");
maxListIter1 = 14; // The maximum index for lyrList for the first iteration

beamSpot = new double[3];
Expand All @@ -231,22 +213,6 @@ public void print() {
// beamSpot[1] = beamPosKal.v[1];
// beamSpot[2] = beamPosKal.v[2];
// }

// Swap axial/stereo in list entries for the top tracker
for (int[] list: lyrList[0]) {
int [] listTop = new int[5];
for (int i=0; i<5; ++i) {
listTop[i] = Swap[list[i]];
}
for (int i=0; i<4; ++i) {
if (listTop[i] > listTop[i+1]) { // Sorting entries. No more than one swap should be necessary.
int tmp = listTop[i];
listTop[i] = listTop[i+1];
listTop[i+1] = tmp;
}
}
lyrList[1].add(listTop);
}
}

public void setLowPhThreshold(double cut) {
Expand Down Expand Up @@ -489,6 +455,10 @@ public void setNumSeedIter1(int num) {
maxListIter1 = n-1;
}

private boolean addStrategy(String strategy) {
return addStrategy(strategy,"top") && addStrategy(strategy,"bottom");
}

// Add a seed search strategy for the bottom or top tracker
public boolean addStrategy(String strategy, String topBottom) {
if (!(topBottom == "top" || topBottom == "bottom")) {
Expand Down
4 changes: 2 additions & 2 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1436,7 +1436,7 @@ private void filterTrack(TrackCandidate tkrCandidate, int lyrBegin, // layer on
int rF;
double [] tRange = {tkrCandidate.tMax - kPar.mxTdif, tkrCandidate.tMin + kPar.mxTdif};
if (prevSite == null) { // For first layer use the initializer state vector
boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made if hitno=-1
boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made on last module of the layer
rF = newSite.makePrediction(sI, null, hitno, tkrCandidate.nTaken <= kPar.mxShared, pickUp, checkBounds, tRange, trial);
if (rF > 0) {
if (m.hits.get(newSite.hitID).tracks.size() > 0) tkrCandidate.nTaken++;
Expand All @@ -1460,7 +1460,7 @@ private void filterTrack(TrackCandidate tkrCandidate, int lyrBegin, // layer on
return;
}
} else {
boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made if hitno=-1
boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made on last module of layer
rF = newSite.makePrediction(prevSite.aF, prevSite.m, hitno, tkrCandidate.nTaken <= kPar.mxShared, pickUp, checkBounds, tRange, trial);
if (rF > 0) {
if (m.hits.get(newSite.hitID).tracks.size() > 0) tkrCandidate.nTaken++;
Expand Down
62 changes: 45 additions & 17 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,9 @@
import hep.aida.IHistogram1D;
import hep.aida.IHistogramFactory;
import hep.physics.vec.Hep3Vector;

/**
* Histograms and plots for Kalman Filter pattern recognition development
* @author Robert Johnson
*
*/
class KalmanPatRecPlots {
private KalmanInterface KI;
Expand Down Expand Up @@ -93,10 +92,9 @@ class KalmanPatRecPlots {

// arguments to histogram1D: name, nbins, min, max
aida.histogram1D("Kalman number of tracks", 10, 0., 10.);
aida.histogram1D("Kalman Track Chi2", 50, 0., 100.);
aida.histogram1D("Kalman Track Chi2, >=12 hits", 50, 0., 100.);
aida.histogram1D("Kalman Track Chi2", 100, 0., 200.);
aida.histogram1D("Kalman Track Chi2, >=12 hits", 100, 0., 200.);
aida.histogram1D("Kalman Track simple Chi2, >=12 hits", 50, 0., 100.);
aida.histogram1D("Kalman Track Chi2 with >=12 good hits", 50, 0., 100.);
aida.histogram2D("number tracks Kalman vs GBL", 20, 0., 5., 20, 0., 5.);
aida.histogram1D("helix chi-squared at origin", 100, 0., 25.);
aida.histogram1D("GBL track chi^2", 50, 0., 100.);
Expand Down Expand Up @@ -184,6 +182,13 @@ class KalmanPatRecPlots {
hnh = aida.histogram1D("MC number hits",15,0.,15.);
hnhf = aida.histogram1D("MC number hits, found",15,0.,15.);
pEff = new Efficiency(40,0.,0.1,"Track efficency vs momentum","momentum (GeV)","efficiency");
aida.histogram1D("Bad/Number of hits on bad tracks", 20, 0., 20.);
aida.histogram1D("Bad/Chi-squared of bad tracks", 100, 0., 200.);
aida.histogram1D("Bad/drho of bad tracks", 50, -8., 8.);
aida.histogram1D("Bad/dz of bad tracks", 50, -4., 4.);
aida.histogram1D("Bad/momentum of bad tracks", 60, 0., 6.);
aida.histogram1D("Bad/Number of MC particles associated", 10, 0., 10.);
aida.histogram1D("Bad/Number of wrong hits on track", 20, 0., 20.);
}

void process(EventHeader event, List<HpsSiSensor> sensors, ArrayList<KalTrack>[] kPatList,
Expand Down Expand Up @@ -419,14 +424,24 @@ void process(EventHeader event, List<HpsSiSensor> sensors, ArrayList<KalTrack>[]
aida.histogram1D("Kalman track time range (ns)").fill(kTk.tMax - kTk.tMin);

// Check the covariance matrix
boolean badCov = false;
Matrix C = new Matrix(kTk.originCovariance());
EigenvalueDecomposition eED= new EigenvalueDecomposition(C);
double [] e = eED.getRealEigenvalues();
for (int i=0; i<5; ++i) {
if (e[i] < 0.) {
logger.warning(String.format("Event %d, eigenvalue %d of covariance is negative!", event.getEventNumber(), i));
System.out.format("Event %d, eigenvalue %d of covariance is negative for track %d!", event.getEventNumber(), i, kTk.ID);
badCov = true;
}
}
if (badCov || kTk.bad) {
aida.histogram1D("Bad/Number of hits on bad tracks").fill(kTk.nHits);
aida.histogram1D("Bad/Chi-squared of bad tracks").fill(kTk.chi2);
aida.histogram1D("Bad/drho of bad tracks").fill(kTk.originHelixParms()[0]);
aida.histogram1D("Bad/dz of bad tracks").fill(kTk.originHelixParms()[3]);
aida.histogram1D("Bad/momentum of bad tracks").fill(pMag);
}
if (kTk.nHits >= 10) {
DMatrixRMaj thisCov = kTk.SiteList.get(0).aS.helix.C;
double ptInv1 = kTk.SiteList.get(0).aS.helix.a.v[2];
Expand Down Expand Up @@ -516,6 +531,7 @@ void process(EventHeader event, List<HpsSiSensor> sensors, ArrayList<KalTrack>[]
}
aida.histogram1D("Kalman track number of shared hits").fill(nShared);
aida.histogram1D("Kalman track number MC particles").fill(mcParts.size());

// Which MC particle is the best match?
int idBest = -1;
int nMatch = 0;
Expand Down Expand Up @@ -548,6 +564,10 @@ void process(EventHeader event, List<HpsSiSensor> sensors, ArrayList<KalTrack>[]
if (!goodHit) nBad++;
}
aida.histogram1D("Kalman number of wrong hits on track").fill(nBad);
if (badCov || kTk.bad) {
aida.histogram1D("Bad/Number of MC particles associated").fill(mcParts.size());
aida.histogram1D("Bad/Number of wrong hits on track").fill(nBad);
}

if (kTk.nHits >= 10) aida.histogram1D("Kalman number of wrong hits on track, >= 10 hits").fill(nBad);
MCParticle mcBest = null;
Expand Down Expand Up @@ -824,21 +844,29 @@ void process(EventHeader event, List<HpsSiSensor> sensors, ArrayList<KalTrack>[]
} //loop on GBL Tracks
} //check if event has GBLTracks

//int [] badEvents = {51753, 52531, 56183, 57958, 58050, 60199, 80324, 83798, 84933, 86351, 88796, 96749, 97230, 102986, 105578, 106654,
// 107191, 108542, 108886, 110453, 120457, 121129, 121311, 121525, 124355, 124910, 127335, 129360, 133951};
int [] badEvents = {};
if (nPlotted < numEvtPlots) { // && sharedHitTrack) {
KI.plotKalmanEvent(outputGnuPlotDir, event, kPatList);
//KI.plotGBLtracks(outputGnuPlotDir, event);
nPlotted++;
}

double maxTruErr = simHitRes(event);
if (maxTruErr < 0.02) {
for (int topBottom=0; topBottom<2; ++topBottom) {
for (KalTrack kTk : kPatList[topBottom]) {
if (kTk.nHits < 12) continue;
aida.histogram1D("Kalman Track Chi2 with >=12 good hits").fill(kTk.chi2);
}
boolean plotIt = false;
if (badEvents.length > 0) {
for (int i=0; i<badEvents.length; ++i) {
if (event.getEventNumber() == badEvents[i]) {
plotIt = true;
break;
}
}
} else {
plotIt = true;
}
if (plotIt) {
KI.plotKalmanEvent(outputGnuPlotDir, event, kPatList);
//KI.plotGBLtracks(outputGnuPlotDir, event);
nPlotted++;
}
}

simHitRes(event);
}

// Make histograms of the MC hit resolution
Expand Down
Loading