Skip to content

Commit ee6e1aa

Browse files
mgignacmgignac
authored
Merge pull request #1159 from JeffersonLab/addisolations
Add isolations to TrackData object for Kalman tracks
2 parents 11e1108 + c0c352d commit ee6e1aa

3 files changed

Lines changed: 90 additions & 43 deletions

File tree

‎tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java‎

Lines changed: 50 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,13 +4,14 @@
44
import java.io.FileNotFoundException;
55
import java.io.PrintWriter;
66
import java.util.ArrayList;
7+
import java.util.List;
78
import java.util.Collections;
89
import java.util.Comparator;
910
import java.util.HashMap;
1011
import java.util.Map;
1112
import java.util.logging.Level;
1213
import java.util.logging.Logger;
13-
14+
import java.util.Comparator;
1415
import org.apache.commons.math.util.FastMath;
1516
import org.ejml.data.DMatrixRMaj;
1617
import org.ejml.dense.row.CommonOps_DDRM;
@@ -514,6 +515,54 @@ public Pair<Double[], Double> unbiasedIntersect(int layer, boolean local) {
514515
}
515516
}
516517

518+
public Pair<Double,Double> getIsoAndT0(int layer){
519+
if (lyrMap == null) {
520+
makeLyrMap();
521+
}
522+
if (lyrMap.containsKey(layer)) {
523+
return getIsoAndT0(lyrMap.get(layer));
524+
}else{
525+
return new Pair<Double, Double>(-999., -999.);
526+
}
527+
}
528+
529+
public Pair<Double, Double> getIsoAndT0(MeasurementSite ms){
530+
double iso=-999.;
531+
double isot0=-999.;
532+
SiModule sensor=ms.m;
533+
List<Measurement> allHits=sensor.hits;
534+
if(allHits.size()>1 && ms.hitID>-1){
535+
Measurement hitOnTrack= allHits.get(ms.hitID);
536+
double hitTime=hitOnTrack.time;
537+
// keep track of sensor position & orientation
538+
// use the sign of the global plane position (vertical = z)
539+
// times the local-->global v-->z element
540+
// to determine if sensor is aligned (+ive v --> away from beam)
541+
// or anti-aligned (+ive v is towards beam)
542+
double vertPos=sensor.p.X().v[2];
543+
double l2gv2z=sensor.R.M[1][2];
544+
int awayFromBeam = (int)Math.signum(vertPos*l2gv2z);
545+
//sort the hits on the module by increasing v
546+
Collections.sort(allHits,Measurement.MeasurementComparatorUp);
547+
// now get the position of nearest hit _away_ from beam
548+
// within 40ns of original hit
549+
int nSteps=1;
550+
int isoID=ms.hitID+awayFromBeam*nSteps;
551+
while(isoID<allHits.size() && isoID>-1){
552+
if(Math.abs(allHits.get(isoID).time-hitTime)<40.0){
553+
iso=Math.abs(allHits.get(isoID).v-hitOnTrack.v);
554+
isot0=allHits.get(isoID).time;
555+
break;
556+
}//otherwise step to the next one
557+
nSteps++;
558+
isoID=ms.hitID+awayFromBeam*nSteps;
559+
}
560+
return new Pair<Double, Double>(iso, isot0);
561+
}else{
562+
return new Pair<Double, Double>(-999., -999.);
563+
}
564+
}
565+
517566
// Returns the unbiased residual for the track at a given layer, together with the variance on that residual
518567
public Pair<Double, Double> unbiasedResidual(MeasurementSite site) {
519568
double resid = -999.;

‎tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java‎

Lines changed: 20 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,6 @@
1717
import org.hps.recon.tracking.TrackData;
1818
import org.hps.recon.tracking.TrackIntersectData;
1919
import org.hps.recon.tracking.TrackResidualsData;
20-
//import org.hps.recon.tracking.KFKinkData;
2120

2221
import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
2322
import org.hps.recon.tracking.MaterialSupervisor.SiStripPlane;
@@ -354,22 +353,16 @@ public void process(EventHeader event) {
354353
List<LCRelation> trackResidualsRelations = new ArrayList<LCRelation>();
355354
List<TrackIntersectData> trackIntersects = new ArrayList<TrackIntersectData>();
356355
List<LCRelation> trackIntersectsRelations = new ArrayList<LCRelation>();
357-
356+
358357
ArrayList<KalTrack>[] kPatList = prepareTrackCollections(event, outputFullTracks, trackDataCollection, trackDataRelations, allClstrs, gblStripClusterDataRelations,trackXKinks,trackXKinksRelations,trackZKinks,trackZKinksRelations,trackResiduals, trackResidualsRelations, trackIntersects, trackIntersectsRelations);
359-
360-
// ArrayList<KalTrack>[] kPatList = prepareTrackCollections(event, outputFullTracks, trackDataCollection, trackDataRelations, allClstrs, gblStripClusterDataRelations, trackResiduals, trackResidualsRelations);
361-
//mg debug why the track data relations (and others I think) are screwed
362-
// for (LCRelation tdRel: trackDataRelations){
363-
// System.out.println(tdRel.getFrom()+" ---> " +tdRel.getTo());
364-
// }
365358

366359
int flag = 1 << LCIOConstants.TRBIT_HITS;
367360
event.put(outputFullTrackCollectionName, outputFullTracks, Track.class, flag);
368361
event.put("KFGBLStripClusterData", allClstrs, GBLStripClusterData.class, flag);
369362
event.put("KFGBLStripClusterDataRelations", gblStripClusterDataRelations, LCRelation.class, flag);
370363
event.put("KFTrackData",trackDataCollection, TrackData.class,0);
371364
event.put("KFTrackDataRelations",trackDataRelations,LCRelation.class,0);
372-
365+
373366
if (addKinks) {
374367
event.put("KFXKink", trackXKinks, TrackResidualsData.class,0);
375368
event.put("KFXKinkRelations", trackXKinksRelations,LCRelation.class,0);
@@ -419,7 +412,8 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
419412
List<TrackResidualsData> trackXKinks, List<LCRelation> trackXKinksRelations,
420413
List<TrackResidualsData> trackZKinks, List<LCRelation> trackZKinksRelations,
421414
List<TrackResidualsData> trackResiduals, List<LCRelation> trackResidualsRelations,
422-
List<TrackIntersectData> trackIntersects, List<LCRelation> trackIntersectsRelations) {
415+
List<TrackIntersectData> trackIntersects, List<LCRelation> trackIntersectsRelations
416+
) {
423417

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

447-
//List<RawTrackerHit> rawhits = event.get(RawTrackerHit.class, "SVTRawTrackerHits");
448-
//if (rawhits == null) {
449-
// logger.log(Level.FINE, String.format("KalmanPatRecDriver.process: the raw hits collection is missing"));
450-
// return null;
451-
//}
452-
453441
int nKalTracks = 0;
454442
for (int topBottom=0; topBottom<2; ++topBottom) {
455443
ArrayList<KalTrack> kPat = kPatList[topBottom];
@@ -472,11 +460,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
472460
Track KalmanTrackHPS = KI.createTrack(kTk, true);
473461
if (KalmanTrackHPS == null) continue;
474462

475-
//pT cut
476-
//double [] hParams_check = kTk.originHelixParms();
477-
//double ptInv_check = hParams_check[2];
478-
//double pt = Math.abs(1./ptInv_check);
479-
480463
outputFullTracks.add(KalmanTrackHPS);
481464

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

503-
//TODO: compute isolations
504-
double qualityArray[] = new double[1];
505-
qualityArray[0] = -1;
506486

507487
//Get the track momentum and convert it into detector frame and float values
508488
Hep3Vector momentum = new BasicHep3Vector(KalmanTrackHPS.getTrackStates().get(0).getMomentum());
@@ -537,13 +517,7 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
537517
//Add the TrackResiduals
538518
List<Integer> layers = new ArrayList<Integer>();
539519
List<Double> residuals = new ArrayList<Double>();
540-
List<Float> sigmas = new ArrayList<Float>();
541-
List<Integer> layersInt = new ArrayList<Integer>();
542-
List<Double> intersect = new ArrayList<Double>();
543-
List<Float> sigmasInt = new ArrayList<Float>();
544-
int uindex = 0;
545-
int vindex = 1;
546-
int windex = 2;
520+
List<Float> sigmas = new ArrayList<Float>();
547521
for (GBLStripClusterData clstr: clstrs) {
548522
Pair<Double,Double> res_and_sigma = kTk.unbiasedResidualMillipede(clstr.getId());
549523
if (res_and_sigma.getSecondElement() > -1.) {
@@ -556,17 +530,30 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
556530
sigmas.add(res_and_sigma.getSecondElement().floatValue());
557531
}
558532
}
533+
//Add the Track Intersections
534+
List<Integer> layersInt = new ArrayList<Integer>();
535+
List<Double> intersect = new ArrayList<Double>();
536+
List<Float> sigmasInt = new ArrayList<Float>();
537+
538+
int uindex = 0;
539+
int vindex = 1;
540+
int windex = 2;
541+
542+
double[] isolationsArray=new double[14];
559543
for(int ilay = 0;ilay<14;ilay++){
560544
Pair<Double[], Double> inter_and_sigma = kTk.unbiasedIntersect(ilay, true);
561545
layersInt.add(ilay);
562546
intersect.add(inter_and_sigma.getFirstElement()[uindex]);
563547
intersect.add(inter_and_sigma.getFirstElement()[vindex]);
564548
intersect.add(inter_and_sigma.getFirstElement()[windex]);
565-
sigmasInt.add(inter_and_sigma.getSecondElement().floatValue());
549+
sigmasInt.add(inter_and_sigma.getSecondElement().floatValue());
550+
//get isolations
551+
Pair<Double,Double> isolation=kTk.getIsoAndT0(ilay);
552+
isolationsArray[ilay]=isolation.getFirstElement();
566553
}//Loop on layers
567554

568555
//Add the Track Data
569-
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY, (float) svtCenter_bFieldY);
556+
TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), isolationsArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY, (float) svtCenter_bFieldY);
570557
trackDataCollection.add(KFtrackData);
571558
trackDataRelations.add(new BaseLCRelation(KFtrackData, KalmanTrackHPS));
572559

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

580-
581567
//Add the Kinks
582568
layers = new ArrayList<Integer>();
583569
List<Double> Xkinks = new ArrayList<Double>();
@@ -598,13 +584,6 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
598584
TrackResidualsData kinkZData = new TrackResidualsData(trackerVolume,layers,Zkinks,sigmas);
599585
trackZKinks.add(kinkZData);
600586
trackZKinksRelations.add(new BaseLCRelation(kinkZData, KalmanTrackHPS));
601-
/*
602-
if (KalmanTrackHPS.getTrackerHits().size() != residuals.size()) {
603-
System.out.println("KalmanPatRecDriver::Residuals consistency check failed.");
604-
System.out.printf("Track has %d hits while I have %d residuals \n", KalmanTrackHPS.getTrackerHits().size(), residuals.size());
605-
}
606-
*/
607-
608587
} // end of loop on tracks
609588
} // end of loop on trackers
610589

‎tracking/src/main/java/org/hps/recon/tracking/kalman/Measurement.java‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
package org.hps.recon.tracking.kalman;
22

33
import java.util.ArrayList;
4-
4+
import java.util.Comparator;
55
/**
66
* Holds a single silicon-strip measurement (single-sided), to interface with the Kalman fit
77
*/
@@ -63,6 +63,10 @@ void addMC(int idx) {
6363
tksMC.add(idx);
6464
}
6565

66+
double getMeasuredValue(){
67+
return v;
68+
};
69+
6670
void print(String s) {
6771
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);
6872
if (tracks.size() == 0) {
@@ -87,4 +91,19 @@ String toString(String s) {
8791
if (rGlobal != null) str = str + rGlobal.toString("global location from MC truth");
8892
return str;
8993
}
94+
95+
// Comparator functions for sorting measurements by measured coordinate value
96+
static Comparator<Measurement> MeasurementComparatorUp = new Comparator<Measurement>() {
97+
public int compare(Measurement s1, Measurement s2) {
98+
double v1 = s1.v;
99+
double v2 = s2.v;
100+
if(v1==v2)
101+
return(0);
102+
else if(v1>v2)
103+
return(1);
104+
else
105+
return(-1);
106+
}
107+
};
108+
90109
}

0 commit comments

Comments
 (0)