Skip to content
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.

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

Iss930 Kalman code for incorporating ECAL energy information into track fits #992

Open
wants to merge 22 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
7ed4f99
Major update of comments; add uniformB option for tracking; added deb…
robertprestonjohnson Nov 18, 2022
51d59ea
New code for ECAL constrained track
robertprestonjohnson Nov 24, 2022
db1c369
documentation update for ECAL constrained tracks
robertprestonjohnson Nov 24, 2022
685b253
Create ReconstructedParticleRefitter.java
Dec 20, 2022
156e361
implement smoothing of tracks that have ECAL energy constraint
robertprestonjohnson Dec 21, 2022
98780df
Merge branch 'iss930' of https://github.com/JeffersonLab/hps-java int…
robertprestonjohnson Dec 21, 2022
cf29750
added hpsHelixIntersect method
robertprestonjohnson Jan 3, 2023
6df8659
added public declaration to the new method hpsHelixIntersect
robertprestonjohnson Jan 3, 2023
30a8fbf
Fixed style violations
Jan 18, 2023
b6fe734
Merge branch 'iss930' of https://github.com/JeffersonLab/hps-java int…
Jan 18, 2023
4391b8b
Fix style check violations
Jan 18, 2023
d676393
recon particle refitting interfacing
robertprestonjohnson Jan 18, 2023
3274444
merge with Normans edits
robertprestonjohnson Jan 18, 2023
2d5f87f
fixed compile and style problems
robertprestonjohnson Jan 19, 2023
935c0a8
debugging
robertprestonjohnson Jan 20, 2023
baf99e3
Refit with energy code complete, although the results are very suspec…
robertprestonjohnson Feb 24, 2023
222b457
Fixed a bug, and now it even seems to work
robertprestonjohnson Mar 3, 2023
0e0f723
refit with and without E constraint seems to work now
robertprestonjohnson May 2, 2023
16befb6
fixed compilation error and removed unused argument in a KalTrack method
robertprestonjohnson May 2, 2023
a93d0ab
update from master
robertprestonjohnson May 3, 2023
c2b5ed1
fixed bugs in hit-removal code
robertprestonjohnson May 11, 2023
3550d95
Merge branch 'master' of https://github.com/JeffersonLab/hps-java int…
robertprestonjohnson May 23, 2023
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
Prev Previous commit
Next Next commit
refit with and without E constraint seems to work now
robertprestonjohnson committed May 2, 2023
commit 0e0f7239ec167572badf5ce33ac292891ad699a1

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
@@ -20,6 +20,7 @@

import org.hps.util.Pair;
import org.lcsim.event.TrackState;
import org.lcsim.geometry.Detector;

/**
* This is for stand-alone testing of the Kalman fit only and is not part of the
@@ -505,7 +506,8 @@ class HelixTest3 { // Program for testing the Kalman fitting code
TkInitial.print("TkInitial: initial helix at the origin");
helixBegin.print("helixBegin: starting helix at layer 1");
RKhelix TkEnd = helixBeginRK;
KalmanInterface KI = new KalmanInterface(kPar, fM);
Detector detect = null;
KalmanInterface KI = new KalmanInterface(kPar, detect, fM);
for (int iTrial = 0; iTrial < nTrials; iTrial++) {
RKhelix Tk = helixBeginRK.copy();
if (verbose) {
@@ -812,7 +814,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code
}
double sigmaE = 0.03 * FastMath.sqrt(Etrue);
double E = Etrue + rnd.nextGaussian() * sigmaE;
KalmanTrack.energyConstraint(E, sigmaE);
KalmanTrack.smoothIt(false, E, sigmaE);

// Check on the covariance matrix
Matrix C = new Matrix(KalmanTrack.originCovariance());
129 changes: 68 additions & 61 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java
Original file line number Diff line number Diff line change
@@ -136,6 +136,12 @@ public class KalTrack {
}
this.SiteList.add(site);
}
if (debug) {
for (MeasurementSite site: SiteList) {
SiModule mod = site.m;
System.out.format("KalTrack: SiModule on layer %d wafer %d, %d hits\n", mod.Layer, mod.detector, mod.hits.size());
}
}

helixAtOrigin = null;
propagated = false;
@@ -587,9 +593,17 @@ public void print(String s) {
for (int i = 0; i < SiteList.size(); i++) {
MeasurementSite site = SiteList.get(i);
SiModule m = site.m;
StateVector aSite = site.aS;
if (aSite == null) aSite = site.aF;
double energy = 0.;
if (aSite != null) {
double tanL = aSite.helix.a.v[4];
double K = aSite.helix.a.v[2];
energy = FastMath.sqrt(1.0+tanL*tanL)/Math.abs(K);
}
int hitID = site.hitID;
System.out.format(" Layer %d, detector %d, stereo=%b, chi^2 inc.=%10.6f ", m.Layer, m.detector, m.isStereo,
site.chi2inc);
System.out.format(" Layer %d, detector %d, stereo=%b, chi^2 inc.=%10.6f E=%8.3f ", m.Layer, m.detector, m.isStereo,
site.chi2inc, energy);
if (hitID >= 0) {
if (site.aS == null) aA = site.aF;
else aA = site.aS;
@@ -1648,8 +1662,7 @@ public int addHits(ArrayList<SiModule> data, double mxResid, double mxChi2inc, d
/**
* Re-fit the track
*
* @param keep true if there might be another recursion after dropping more
* hits
* @param keep true if there might be another recursion after dropping more hits
* @return true if the refit was successful
*/
public boolean fit(boolean keep) {
@@ -1787,10 +1800,11 @@ public boolean fit(boolean keep) {
* momentum measurement). Then smooth back to the first layer and propagate
* to the origin.
*
* @param E ECAL energy
* @param sigmaE ECAL energy uncertainty
* @param eConstraint true to include the energy constraint
* @param E ECAL energy
* @param sigmaE ECAL energy uncertainty
*/
void energyConstraint(double E, double sigmaE) {
void smoothIt(boolean eConstraint, double E, double sigmaE) {
// The prediction step from the last tracker site to the ECAL has F=1,
// since we don't alter an helix parameters in the prediction.
// The HelixState.energyConstrained method then uses the Kalman weighted
@@ -1800,71 +1814,64 @@ void energyConstraint(double E, double sigmaE) {
// site is exactly what is returned by the energyConstrained method. The
// smoothing back to the first tracker site can then proceed as usual.
//boolean debug = true;
final boolean noConstraint = false; // for debugging purposes, should just do a normal Kalman smoothing
if (debug) System.out.format("Entering KalTrack.smoothIt, eConstraint=%b, E=%f, sigmaE=%f\n", eConstraint, E, sigmaE);
int idxLast = SiteList.size() - 1;
MeasurementSite lastSite = SiteList.get(idxLast);
HelixState energyConstrainedHelix = lastSite.aF.helix.energyConstrained(E, sigmaE);
if (debug) System.out.format("KalTrack:energyConstraint kappa=%9.4f, kappaE=%9.4f\n",lastSite.aF.helix.a.v[2],energyConstrainedHelix.a.v[2]);
lastSite.energyConstrained = !noConstraint;
lastSite.aS = new StateVector(lastSite.aF.kLow, lastSite.aF.uniformB); // Blank new state vector
if (noConstraint) {
lastSite.aS.helix = lastSite.aF.helix;
lastSite.energyConstrained = eConstraint;
if (!eConstraint) {
lastSite.aS = lastSite.aF;
this.chi2 = lastSite.chi2inc;
} else {
if (debug) {
double kappa = lastSite.aF.helix.a.v[2];
double tanl = lastSite.aF.helix.a.v[4];
double ePredict = FastMath.sqrt(1.0 + tanl * tanl) / Math.abs(kappa);
System.out.format("KalTrack.smoothIt at last layer, E=%9.4f\n", ePredict);
}
lastSite.aS = new StateVector(lastSite.aF.kLow, lastSite.aF.uniformB); // Blank new state vector
HelixState energyConstrainedHelix = lastSite.aF.helix.energyConstrained(E, sigmaE);
lastSite.aS.helix = energyConstrainedHelix;
lastSite.aS.kUp = lastSite.aF.kUp;
lastSite.aS.F = lastSite.aF.F; // Don't deep copy the F matrix
lastSite.aS.mPred = lastSite.aF.mPred;
lastSite.aS.R = lastSite.aF.R;
lastSite.aS.r = lastSite.aF.r;
lastSite.aS.K = lastSite.aF.K;
lastSite.smoothed = true;
lastSite.chi2inc = (lastSite.aS.r * lastSite.aS.r) / lastSite.aS.R;

// Get the residual of the prediction at the ECAL
double kappa = energyConstrainedHelix.a.v[2];
double tanl = energyConstrainedHelix.a.v[4];
double ePredict = FastMath.sqrt(1.0 + tanl * tanl) / Math.abs(kappa);
double chi = (E - ePredict) / sigmaE;
this.chi2_Econstraint = lastSite.chi2inc + chi * chi;
if (debug) {
System.out.format("KalTrack:smoothIt E=%9.4f, Epredict=%9.4f, sigmaE=%9.4f, chi=%9.4f\n",
E,ePredict,sigmaE,chi);
System.out.format(" constrained helix = %s\n", energyConstrainedHelix.a.toString());
}
}
lastSite.aS.kUp = lastSite.aF.kUp;
lastSite.aS.F = lastSite.aF.F; // Don't deep copy the F matrix
lastSite.aS.mPred = lastSite.aF.mPred;
lastSite.aS.R = lastSite.aF.R;
lastSite.aS.r = lastSite.aF.r;
lastSite.aS.K = lastSite.aF.K;
lastSite.chi2incE = (lastSite.aS.r * lastSite.aS.r) / lastSite.aS.R;
// Get the residual of the prediction at the ECAL
double kappa = energyConstrainedHelix.a.v[2];
double tanl = energyConstrainedHelix.a.v[4];
double ePredict = FastMath.sqrt(1.0 + tanl * tanl) / Math.abs(kappa);
double chi = (E - ePredict) / sigmaE;
if (debug) System.out.format("KalTrack:energyConstraint E=%9.4f, Epredict=%9.4f, sigmaE=%9.4f, chi=%9.4f\n",E,ePredict,sigmaE,chi);
if (noConstraint) {
this.chi2_Econstraint = lastSite.chi2incE;
} else {
this.chi2_Econstraint = lastSite.chi2incE + chi * chi;
}
lastSite.smoothed = true;
MeasurementSite nS = lastSite;
for (int idx = idxLast - 1; idx >= 0; --idx) {
MeasurementSite thisSite = SiteList.get(idx);
thisSite.aS = thisSite.aF.smooth(nS.aS, nS.aP);
if (thisSite.hitID < 0) {
thisSite.energyConstrained = !noConstraint;
continue;
}
Measurement hit = thisSite.m.hits.get(thisSite.hitID);
double V = hit.sigma * hit.sigma;
double phiS = thisSite.aS.helix.planeIntersect(thisSite.m.p);

if (Double.isNaN(phiS)) { // This should almost never happen!
logger.log(Level.FINE, "KalTrack.energyConstraint: no intersection of helix with the plane exists.");
continue;
}
thisSite.aS.mPred = thisSite.h(thisSite.aS, thisSite.m, phiS);
thisSite.aS.r = hit.v - thisSite.aS.mPred;
if (tempV == null) {
tempV = new DMatrixRMaj(5, 1);
thisSite.smooth(nS);
thisSite.energyConstrained = eConstraint;
thisSite.smoothed = true;
if (eConstraint) {
this.chi2_Econstraint += thisSite.chi2inc;
} else {
this.chi2 += thisSite.chi2inc;
}
CommonOps_DDRM.mult(thisSite.aS.helix.C, thisSite.H, tempV);
thisSite.aS.R = V - CommonOps_DDRM.dot(thisSite.H, tempV);
if (thisSite.aS.R < 0) {
if (KalmanPatRecHPS.negativeCov(thisSite.aS.helix.C)) {
if (debug) {
System.out.format("KalTrack.energyConstraint, measurement covariance %12.4e is negative\n", thisSite.aS.R);
System.out.format("KalTrack: event %d, ID %d, negative covariance after smoothing at layer %d\n",
eventNumber, ID, thisSite.m.Layer);
}
//aS.print("the smoothed state");
//nS.print("the next site in the chain");
thisSite.aS.R = 0.25 * V; // A negative covariance makes no sense, hence this fudge
bad = true;
KalmanPatRecHPS.fixCov(thisSite.aS.helix.C, thisSite.aS.helix.a);
}

thisSite.chi2incE = (thisSite.aS.r * thisSite.aS.r) / thisSite.aS.R;
this.chi2_Econstraint += thisSite.chi2incE;
thisSite.energyConstrained = !noConstraint;
nS = thisSite;
}
Vec beamSpot = new Vec(3, kPar.beamSpot);
@@ -1874,7 +1881,7 @@ void energyConstraint(double E, double sigmaE) {
Plane originPlane = new Plane(beamSpot, new Vec(0., 1., 0.));
arcLengthE = new double[1];
helixAtOrigin = SiteList.get(0).aS.helix.propagateRungeKutta(originPlane, yScat, XLscat, SiteList.get(0).m.Bfield, arcLengthE);
energyConstrained = true;
energyConstrained = eConstraint;
}

/**
Original file line number Diff line number Diff line change
@@ -174,7 +174,7 @@ public void detectorChanged(Detector det) {

KalmanParams kPar = new KalmanParams();
kPar.setUniformB(uniformB);
KI = new KalmanInterface(kPar, fm);
KI = new KalmanInterface(kPar, det, fm);
KI.createSiModules(detPlanes);

System.out.format("KalmanDriver: the B field is assumed uniform? %b\n", uniformB);

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
@@ -69,7 +69,7 @@ public void detectorChanged(Detector det) {
KalmanParams kPar = new KalmanParams();
kPar.print();

KI = new KalmanInterface(kPar, fm);
KI = new KalmanInterface(kPar, det, fm);
KI.createSiModules(detPlanes);

}
Original file line number Diff line number Diff line change
@@ -173,7 +173,7 @@ public KalmanParams() {
lowPhThresh = 0.25; // Residual improvement ratio necessary to use a low-ph hit instead of high-ph
seedCompThr = 0.05; // Remove SeedTracks with all Helix params within relative seedCompThr . If -1 do not apply duplicate removal
eRes = new double[2];
eRes[0] = 10.0; // Cal energy resolution parameters in % sigmaE = eRes[0]/sqrt(E) + eRes[1]
eRes[0] = 5.0; // Cal energy resolution parameters in % sigmaE = eRes[0]/sqrt(E) + eRes[1]
eRes[1] = 1.0;

// Load the default search strategies
@@ -239,6 +239,13 @@ public void setEnergyRes(double a, double b) {
if (b > 0.) eRes[1] = b;
logger.config(String.format("Setting CAL energy resolution to %8.2f/sqrt(E) + %8.2f", eRes[0], eRes[1]));
}
public double getEres(int i) {
if (i<0 || i>1) {
logger.warning(String.format("KalmanParams: invalid eRes index %d\n", i));
return eRes[0];
}
return eRes[i];
}

public void setUniformB(boolean input) {
logger.config(String.format("Setting the field to be uniform? %b", input));
Original file line number Diff line number Diff line change
@@ -374,7 +374,7 @@ public void detectorChanged(Detector det) {
logger.config("KalmanPatRecDriver: done with configuration changes.");
kPar.print();

KI = new KalmanInterface(kPar, fm);
KI = new KalmanInterface(kPar, det, fm);
KI.setSiHitsLimit(siHitsLimit);
KI.createSiModules(detPlanes);
decoder = det.getSubdetector("Tracker").getIDDecoder();
@@ -456,10 +456,8 @@ public int compare(TrackerHit o1, TrackerHit o2) {
* @param event input the header for this event
* @param outputFullTracks output the list of HPS tracks
* @param trackDataCollection output list of data that go with the tracks
* @param trackDataRelations output the relations between tracks and the
* data
* @param allClstrs output all the clusters needed for refitting the Kalman
* tracks by GBL
* @param trackDataRelations output the relations between tracks and the data
* @param allClstrs output all the clusters needed for refitting the Kalman tracks by GBL
* @param gblStripClusterDataRelations output relations for the clusters
* @param trackResiduals output the residuals for hits on Kalman tracks
* @param trackResidualsRelations output relations for the residuals
@@ -476,7 +474,7 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
}

long startTime = System.nanoTime();
ArrayList<KalTrack>[] kPatList = KI.KalmanPatRec(event, decoder);
ArrayList<KalTrack>[] kPatList = KI.KalmanPatRec(event);
long endTime = System.nanoTime();
double runTime = (double) (endTime - startTime) / 1000000.;
executionTime += runTime;
Original file line number Diff line number Diff line change
@@ -23,7 +23,6 @@ class MeasurementSite {
boolean smoothed; // True if the smoothed state vector has been built
boolean energyConstrained; // True if the smoothed state vector is energy constrained
double chi2inc; // chi^2 increment for this site
double chi2incE; // chi^2 increment for the energy-constrainted track
DMatrixRMaj H; // Derivatives of the transformation from state vector to measurement
double arcLength; // Arc length from the previous measurement
private double conFac; // Conversion from B to alpha
@@ -76,7 +75,7 @@ String toString(String s) {
Vec tB = Bfield.unitVec();
str=str+String.format(" Magnetic field strength=%10.6f; alpha=%10.6f\n", B, alpha);
str = str + tB.toString("magnetic field direction") + "\n";
str=str+String.format(" chi^2 increment=%12.4e; E constrained=%12.4E\n", chi2inc, chi2incE);
str=str+String.format(" chi^2 increment=%12.4e\n", chi2inc);
str=str+String.format(" x scattering angle=%10.8f, y scattering angle=%10.8f\n", scatX(), scatZ());
if (predicted) str = str + aP.toString("predicted");
if (filtered) str = str + aF.toString("filtered");
@@ -121,7 +120,6 @@ String toString(String s) {
double sp = 0.002; // Estar collision stopping power for electrons in silicon at about a GeV, in GeV cm2/g
dEdx = -0.1 * sp * rho; // in GeV/mm
chi2inc = 0.;
chi2incE = 0.;
H = new DMatrixRMaj(5,1);
if (!initialized) {
tempV = new DMatrixRMaj(5,1);
@@ -557,7 +555,6 @@ boolean removeHit() {
if (hitID < 0) return false;
hitID = -1;
chi2inc = 0.;
chi2incE = 0.;
smoothed = false;
filtered = false;
return true;
Original file line number Diff line number Diff line change
@@ -12,7 +12,7 @@
class SiModule {
int Layer; // Tracker layer number, or a negative integer for a dummy layer added just for stepping in a
// non-uniform field
int detector; // Detector number within the layer
int detector; // Detector or module number within the layer
int millipedeID; // ID used by millipede for alignment
ArrayList<Measurement> hits; // Hits ordered by coordinate value, from minimum to maximum
Plane p; // Orientation and offset of the detector measurement plane in global coordinates