Skip to content

Add isolations to TrackData object for Kalman tracks #1159

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
mgignac merged 7 commits into from
Dec 8, 2025
Merged

Add isolations to TrackData object for Kalman tracks #1159

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
36f0523
add isolations as a generic object and save to event
bloodyyugo Dec 3, 2025
400cdc4
use the double array in TrackData to store the isolations instead
bloodyyugo Dec 3, 2025
cc4bd1f
clean up commented code
bloodyyugo Dec 3, 2025
1e3a722
Update condition to check allHits size in KalTrack
bloodyyugo Dec 5, 2025
4967fff
Remove unused variable hpsSensorPos
bloodyyugo Dec 5, 2025
3dafcbc
require the iso hit be within 40ns of hit on track
bloodyyugo Dec 5, 2025
c0c352d
merge in changes
bloodyyugo Dec 5, 2025
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
51 changes: 50 additions & 1 deletion 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 @@ -4,13 +4,14 @@
import java.io.FileNotFoundException;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.List;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.Map;
import java.util.logging.Level;
import java.util.logging.Logger;

import java.util.Comparator;
import org.apache.commons.math.util.FastMath;
import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
Expand Down Expand Up @@ -514,6 +515,54 @@ public Pair<Double[], Double> unbiasedIntersect(int layer, boolean local) {
}
}

public Pair<Double,Double> getIsoAndT0(int layer){
if (lyrMap == null) {
makeLyrMap();
}
if (lyrMap.containsKey(layer)) {
return getIsoAndT0(lyrMap.get(layer));
}else{
return new Pair<Double, Double>(-999., -999.);
}
}

public Pair<Double, Double> getIsoAndT0(MeasurementSite ms){
double iso=-999.;
double isot0=-999.;
SiModule sensor=ms.m;
List<Measurement> allHits=sensor.hits;
if(allHits.size()>1 && ms.hitID>-1){
Measurement hitOnTrack= allHits.get(ms.hitID);
double hitTime=hitOnTrack.time;
// keep track of sensor position & orientation
// use the sign of the global plane position (vertical = z)
// times the local-->global v-->z element
// to determine if sensor is aligned (+ive v --> away from beam)
// or anti-aligned (+ive v is towards beam)
Copy link
Collaborator

@sarahgaiser sarahgaiser Dec 5, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't fully understand this comment. What do you mean by local-->global v-->z element and what is +ive v

Copy link
Contributor Author

@bloodyyugo bloodyyugo Dec 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The local-->global refers to the rotation matrix R (local i.e. sensor frame to the "global" frame for kalman tracking) and v-->z element is the element in that matrix (the z component of the v direction). +ive == positive.

double vertPos=sensor.p.X().v[2];
double l2gv2z=sensor.R.M[1][2];
int awayFromBeam = (int)Math.signum(vertPos*l2gv2z);
//sort the hits on the module by increasing v
Collections.sort(allHits,Measurement.MeasurementComparatorUp);
// now get the position of nearest hit _away_ from beam
// within 40ns of original hit
int nSteps=1;
int isoID=ms.hitID+awayFromBeam*nSteps;
while(isoID<allHits.size() && isoID>-1){
if(Math.abs(allHits.get(isoID).time-hitTime)<40.0){
iso=Math.abs(allHits.get(isoID).v-hitOnTrack.v);
isot0=allHits.get(isoID).time;
break;
}//otherwise step to the next one
nSteps++;
isoID=ms.hitID+awayFromBeam*nSteps;
}
return new Pair<Double, Double>(iso, isot0);
}else{
return new Pair<Double, Double>(-999., -999.);
}
}

// Returns the unbiased residual for the track at a given layer, together with the variance on that residual
public Pair<Double, Double> unbiasedResidual(MeasurementSite site) {
double resid = -999.;
Expand Down
61 changes: 20 additions & 41 deletions tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
import org.hps.recon.tracking.TrackData;
import org.hps.recon.tracking.TrackIntersectData;
import org.hps.recon.tracking.TrackResidualsData;
//import org.hps.recon.tracking.KFKinkData;

import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
import org.hps.recon.tracking.MaterialSupervisor.SiStripPlane;
Expand Down Expand Up @@ -354,22 +353,16 @@ public void process(EventHeader event) {
List<LCRelation> trackResidualsRelations = new ArrayList<LCRelation>();
List<TrackIntersectData> trackIntersects = new ArrayList<TrackIntersectData>();
List<LCRelation> trackIntersectsRelations = new ArrayList<LCRelation>();

ArrayList<KalTrack>[] kPatList = prepareTrackCollections(event, outputFullTracks, trackDataCollection, trackDataRelations, allClstrs, gblStripClusterDataRelations,trackXKinks,trackXKinksRelations,trackZKinks,trackZKinksRelations,trackResiduals, trackResidualsRelations, trackIntersects, trackIntersectsRelations);

// ArrayList<KalTrack>[] kPatList = prepareTrackCollections(event, outputFullTracks, trackDataCollection, trackDataRelations, allClstrs, gblStripClusterDataRelations, trackResiduals, trackResidualsRelations);
//mg debug why the track data relations (and others I think) are screwed
// for (LCRelation tdRel: trackDataRelations){
// System.out.println(tdRel.getFrom()+" ---> " +tdRel.getTo());
// }

int flag = 1 << LCIOConstants.TRBIT_HITS;
event.put(outputFullTrackCollectionName, outputFullTracks, Track.class, flag);
event.put("KFGBLStripClusterData", allClstrs, GBLStripClusterData.class, flag);
event.put("KFGBLStripClusterDataRelations", gblStripClusterDataRelations, LCRelation.class, flag);
event.put("KFTrackData",trackDataCollection, TrackData.class,0);
event.put("KFTrackDataRelations",trackDataRelations,LCRelation.class,0);

if (addKinks) {
event.put("KFXKink", trackXKinks, TrackResidualsData.class,0);
event.put("KFXKinkRelations", trackXKinksRelations,LCRelation.class,0);
Expand Down Expand Up @@ -419,7 +412,8 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
List<TrackResidualsData> trackXKinks, List<LCRelation> trackXKinksRelations,
List<TrackResidualsData> trackZKinks, List<LCRelation> trackZKinksRelations,
List<TrackResidualsData> trackResiduals, List<LCRelation> trackResidualsRelations,
List<TrackIntersectData> trackIntersects, List<LCRelation> trackIntersectsRelations) {
List<TrackIntersectData> trackIntersects, List<LCRelation> trackIntersectsRelations
) {

int evtNumb = event.getEventNumber();
String stripHitInputCollectionName = "StripClusterer_SiTrackerHitStrip1D";
Expand All @@ -444,12 +438,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
nEvents++;
logger.log(Level.FINE,"KalmanPatRecDriver.process: run time for pattern recognition at event " + evtNumb + " is " + runTime + " milliseconds");

//List<RawTrackerHit> rawhits = event.get(RawTrackerHit.class, "SVTRawTrackerHits");
//if (rawhits == null) {
// logger.log(Level.FINE, String.format("KalmanPatRecDriver.process: the raw hits collection is missing"));
// return null;
//}

int nKalTracks = 0;
for (int topBottom=0; topBottom<2; ++topBottom) {
ArrayList<KalTrack> kPat = kPatList[topBottom];
Expand All @@ -472,11 +460,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
Track KalmanTrackHPS = KI.createTrack(kTk, true);
if (KalmanTrackHPS == null) continue;

//pT cut
//double [] hParams_check = kTk.originHelixParms();
//double ptInv_check = hParams_check[2];
//double pt = Math.abs(1./ptInv_check);

outputFullTracks.add(KalmanTrackHPS);

List<GBLStripClusterData> clstrs = KI.createGBLStripClusterData(kTk);
Expand All @@ -500,9 +483,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
//if tanLamda<0 set bottom
if (KalmanTrackHPS.getTrackStates().get(0).getTanLambda() < 0) trackerVolume = 1;

//TODO: compute isolations
double qualityArray[] = new double[1];
qualityArray[0] = -1;

//Get the track momentum and convert it into detector frame and float values
Hep3Vector momentum = new BasicHep3Vector(KalmanTrackHPS.getTrackStates().get(0).getMomentum());
Expand Down Expand Up @@ -537,13 +517,7 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
//Add the TrackResiduals
List<Integer> layers = new ArrayList<Integer>();
List<Double> residuals = new ArrayList<Double>();
List<Float> sigmas = new ArrayList<Float>();
List<Integer> layersInt = new ArrayList<Integer>();
Copy link
Collaborator

@sarahgaiser sarahgaiser Dec 5, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason why this chunk of code is moved down a few lines? It is hard to parse what is going on with the different loops and if statements.

Copy link
Contributor Author

@bloodyyugo bloodyyugo Dec 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The first loop is of hits on track (so, only layers with a HOT are included) and the second one is over all layers. I separated the intersections list instantiations down to just before the loop over all layers because that's where it's filled.

List<Double> intersect = new ArrayList<Double>();
List<Float> sigmasInt = new ArrayList<Float>();
int uindex = 0;
int vindex = 1;
int windex = 2;
List<Float> sigmas = new ArrayList<Float>();
for (GBLStripClusterData clstr: clstrs) {
Pair<Double,Double> res_and_sigma = kTk.unbiasedResidualMillipede(clstr.getId());
if (res_and_sigma.getSecondElement() > -1.) {
Expand All @@ -556,17 +530,30 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
sigmas.add(res_and_sigma.getSecondElement().floatValue());
}
}
//Add the Track Intersections
List<Integer> layersInt = new ArrayList<Integer>();
List<Double> intersect = new ArrayList<Double>();
List<Float> sigmasInt = new ArrayList<Float>();

int uindex = 0;
int vindex = 1;
int windex = 2;

double[] isolationsArray=new double[14];
for(int ilay = 0;ilay<14;ilay++){
Pair<Double[], Double> inter_and_sigma = kTk.unbiasedIntersect(ilay, true);
layersInt.add(ilay);
intersect.add(inter_and_sigma.getFirstElement()[uindex]);
intersect.add(inter_and_sigma.getFirstElement()[vindex]);
intersect.add(inter_and_sigma.getFirstElement()[windex]);
sigmasInt.add(inter_and_sigma.getSecondElement().floatValue());
sigmasInt.add(inter_and_sigma.getSecondElement().floatValue());
//get isolations
Pair<Double,Double> isolation=kTk.getIsoAndT0(ilay);
isolationsArray[ilay]=isolation.getFirstElement();
}//Loop on layers

//Add the Track Data
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY, (float) svtCenter_bFieldY);
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), isolationsArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY, (float) svtCenter_bFieldY);
trackDataCollection.add(KFtrackData);
trackDataRelations.add(new BaseLCRelation(KFtrackData, KalmanTrackHPS));

Expand All @@ -577,7 +564,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
trackIntersects.add(intersectData);
trackIntersectsRelations.add(new BaseLCRelation(intersectData, KalmanTrackHPS));


//Add the Kinks
layers = new ArrayList<Integer>();
List<Double> Xkinks = new ArrayList<Double>();
Expand All @@ -598,13 +584,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
TrackResidualsData kinkZData = new TrackResidualsData(trackerVolume,layers,Zkinks,sigmas);
trackZKinks.add(kinkZData);
trackZKinksRelations.add(new BaseLCRelation(kinkZData, KalmanTrackHPS));
/*
if (KalmanTrackHPS.getTrackerHits().size() != residuals.size()) {
System.out.println("KalmanPatRecDriver::Residuals consistency check failed.");
System.out.printf("Track has %d hits while I have %d residuals \n", KalmanTrackHPS.getTrackerHits().size(), residuals.size());
}
*/

} // end of loop on tracks
} // end of loop on trackers

Expand Down
21 changes: 20 additions & 1 deletion tracking/src/main/java/org/hps/recon/tracking/kalman/Measurement.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
package org.hps.recon.tracking.kalman;

import java.util.ArrayList;

import java.util.Comparator;
/**
* Holds a single silicon-strip measurement (single-sided), to interface with the Kalman fit
*/
Expand Down Expand Up @@ -63,6 +63,10 @@ void addMC(int idx) {
tksMC.add(idx);
}

double getMeasuredValue(){
return v;
};

void print(String s) {
System.out.format("Measurement %s: Measurement value=%10.5f+-%8.6f; xStrip=%7.2f, MC truth=%10.5f; t=%8.3f; E=%8.3f", s, v, sigma, x, vTrue, time, energy);
if (tracks.size() == 0) {
Expand All @@ -87,4 +91,19 @@ String toString(String s) {
if (rGlobal != null) str = str + rGlobal.toString("global location from MC truth");
return str;
}

// Comparator functions for sorting measurements by measured coordinate value
static Comparator<Measurement> MeasurementComparatorUp = new Comparator<Measurement>() {
public int compare(Measurement s1, Measurement s2) {
double v1 = s1.v;
double v2 = s2.v;
if(v1==v2)
return(0);
else if(v1>v2)
return(1);
else
return(-1);
}
};

}
Loading