Skip to content

Commit c733f6d

Browse files
bloodyyugobloodyyugo
committed
add hole accounting to subdetectorHitNumbers as -1
1 parent 5686a05 commit c733f6d

3 files changed

Lines changed: 97 additions & 26 deletions

File tree

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

Lines changed: 29 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -295,27 +295,49 @@ public Pair<Double[], Double> unbiasedIntersect(MeasurementSite site, boolean lo
295295
Vec intPnt = HelixState.atPhi(site.aS.helix.X0, aStar, phiInt, site.aS.helix.alpha);
296296
// Transform the intersection back into the global coordinates
297297
Vec globalInt = site.aS.helix.toGlobal(intPnt);
298+
// System.out.println("KalTrack::unbiasedIntersect global x = "+globalInt.v[0]+
299+
// " y = "+globalInt.v[1]+
300+
// " z = "+globalInt.v[2]);
301+
298302
if (!local) {
299-
new Pair<>(new Double[]{globalInt.v[0], globalInt.v[1], globalInt.v[2]}, 999); //I didn't include global variance since I'm lazy
303+
new Pair<>(new Double[]{globalInt.v[0], globalInt.v[1], globalInt.v[2]}, 999); //I didn't include global variance since I'm lazy
300304
} // Transform the intersection point to the local sensor system
301305
Vec localInt = site.m.toLocal(globalInt);
302306
CommonOps_DDRM.mult(Cstar, site.H, tempV);
303307
Double varUmeas = CommonOps_DDRM.dot(site.H, tempV);
304308
return new Pair<>(new Double[]{localInt.v[0], localInt.v[1], localInt.v[2]}, varUmeas);
305309
} else {
306-
if(debug)System.out.println("no phi-interect with this layer found!");
310+
if(debug)System.out.println("KalTrack::unbiasedIntersect no phi-interect with this layer found!");
307311
}
308312

309313
return new Pair<>(new Double[]{666., 666., 666.}, 666.);
310314
}
311315

316+
public boolean isTrackHole(int layer){
317+
if (lyrMap == null) {
318+
makeLyrMap();
319+
}
320+
if (lyrMap.containsKey(layer)) {
321+
MeasurementSite ms=lyrMap.get(layer);
322+
return isTrackHole(ms);
323+
}else{
324+
return false;
325+
}
326+
}
327+
public boolean isTrackHole(MeasurementSite ms){
328+
if(ms.hitID >=0 && ms.m.hits.size()>0)
329+
return false; //track has a hit in this layer
330+
Pair<Double[], Double> lInt=unbiasedIntersect(ms,true);
331+
boolean trkInSensor=ms.isInSensor(lInt.getFirstElement());
332+
return trkInSensor;
333+
}
312334
private void makeLyrMap() {
313335
lyrMap = new HashMap<Integer, MeasurementSite>(nHits);
314336
for (MeasurementSite site : SiteList) {
315337
lyrMap.put(site.m.Layer, site);
316338
}
317339
}
318-
340+
319341
private void makeMillipedeMap() {
320342
millipedeMap = new HashMap<Integer, MeasurementSite>(nHits);
321343
for (MeasurementSite site : SiteList) {
@@ -512,7 +534,9 @@ public Pair<Double, Double> unbiasedResidual(MeasurementSite site) {
512534

513535
CommonOps_DDRM.mult(Cstar, site.H, tempV);
514536
varResid = sigma * sigma + CommonOps_DDRM.dot(site.H, tempV);
515-
}
537+
}else{
538+
System.out.println("KalTrack::unbiasedResidual phiInt is NaN");
539+
}
516540
}
517541
return new Pair<Double, Double>(resid, varResid);
518542
}
@@ -908,26 +932,7 @@ public HelixState originConstraint(double[] vtx, double[][] vtxCov) {
908932
Vec v = helixAtOrigin.toLocal(new Vec(3, vtx));
909933
SquareMatrix Cov = helixAtOrigin.Rot.rotate(new SquareMatrix(3, vtxCov));
910934
Vec X0 = helixAtOrigin.X0;
911-
double phi = phiDOCA(helixAtOrigin.a, v, X0, alpha);
912-
/* if (debug) { // Test the DOCA algorithm
913-
Vec rDoca = HelixState.atPhi(X0, helixAtOrigin.a, phi, alpha);
914-
System.out.format("originConstraint: phi of DOCA=%10.5e\n", phi);
915-
rDoca.print(" DOCA point");
916-
double doca = rDoca.dif(v).mag();
917-
System.out.format(" Minimum doca=%10.7f\n", doca);
918-
for (double p=phi-0.0001; p<phi+0.0001; p += 0.000001) {
919-
rDoca = HelixState.atPhi(X0, helixAtOrigin.a, p, alpha);
920-
doca = rDoca.dif(v).mag();
921-
System.out.format(" phi=%10.5e, doca=%10.7f\n", p,doca);
922-
}
923-
924-
double delPhi = 0.00001;
925-
double f = fDOCA(phi, helixAtOrigin.a, v, X0, alpha);
926-
double df1 = fDOCA(phi + delPhi, helixAtOrigin.a, v, X0, alpha) - f;
927-
double deriv = dfDOCAdPhi(phi, helixAtOrigin.a, v, X0, alpha);
928-
double df2 = deriv * delPhi;
929-
System.out.format("Test of fDOCA derivative: df exact = %11.7f; df from derivative = %11.7f\n", df1, df2);
930-
}*/
935+
double phi = phiDOCA(helixAtOrigin.a, v, X0, alpha);
931936
double[][] H = buildH(helixAtOrigin.a, v, X0, phi, alpha);
932937
Vec pntDOCA = HelixState.atPhi(X0, helixAtOrigin.a, phi, alpha);
933938
if (debug) {

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

Lines changed: 60 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,7 @@ public class KalmanInterface {
101101
private int nBigEvents;
102102
private int eventNumber;
103103
private static double target_pos = -999.9;
104+
private static boolean addTrackStateAtTarget = false;
104105
private double[] beamPosition = null;
105106

106107
private static final boolean debug = false;
@@ -651,11 +652,32 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
651652
MeasurementSite site = kT.SiteList.get(i);
652653
ts = null;
653654
int loc = TrackState.AtOther;
655+
//add holes to hitpattern
656+
boolean isTrkHole=kT.isTrackHole(site);
657+
if(isTrkHole){
658+
Array.set(hitPattern,site.m.Layer,-1);
659+
}
660+
//HpsSiSensor hssd = (HpsSiSensor) moduleMap.get(site.m).getSensor();
661+
//int lay = hssd.getMillepedeId();
662+
// System.out.printf("ssp id %d \n", hssd.getMillepedeId());
654663

655664
if (i == firstHit_idx) {
656665
loc = TrackState.AtFirstHit;
657666
} else if (i == lastHit_idx)
658-
loc = TrackState.AtLastHit;
667+
loc = TrackState.AtLastHit;
668+
669+
/*
670+
//DO Not att the missing layer track states yet.
671+
if (storeTrackStates) {
672+
for (int k = 1; k < lay - prevID; k++) {
673+
// uses new lcsim constructor
674+
BaseTrackState dummy = new BaseTrackState(dummyCounter);
675+
newTrack.getTrackStates().add(dummy);
676+
dummyCounter--;
677+
}
678+
prevID = lay;
679+
}
680+
*/
659681

660682
if (loc == TrackState.AtFirstHit || loc == TrackState.AtLastHit || storeTrackStates) {
661683
ts = createTrackState(site, loc, true, saveTrackStateAtIntercept);
@@ -676,14 +698,40 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
676698
double BAtEcal = BfieldAtEcal.mag();
677699
double alphaAtEcal = conFac/ BAtEcal;
678700

701+
/*
702+
703+
DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal));
704+
double[] ecalParams = helixAtEcal.a.v.clone();
705+
double[] ecalLCSimParams = getLCSimParams(ecalParams, alphaAtEcal);
706+
double[] ecalLCSimCov = getLCSimCov(ecalCov, alphaAtEcal).asPackedArray(true);
707+
double[] refAtEcal = localKalToHps(helixAtEcal.origin);
708+
if(debug)System.out.println(this.getClass().getName()+":: reference to TrackState @ ecal = "+
709+
+refAtEcal[0]+", "+refAtEcal[1]+","+ refAtEcal[2]);
710+
TrackState ts_ecal_helix = new BaseTrackState(ecalLCSimParams,refAtEcal, ecalLCSimCov, TrackState.AtCalorimeter, BAtEcal);
711+
if(debug)System.out.println(this.getClass().getName()+":: Uncorrected track state: curvature = "+ts_ecal_helix.getOmega()
712+
+" bField = "+ts_ecal_helix.getBLocal()+" momentum z = "+ts_ecal_helix.getMomentum()[0]);
713+
System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions ");
714+
System.out.println(ts_ecal_helix.toString());
715+
*/
716+
//newTrack.getTrackStates().add(ts_ecal_helix);
717+
718+
679719
// this is how createTrackState (above) goes from measurement to TrackState
680720
// from the measurements. It calls toHPShelix.
681721
// the helix state here must already be propagated
682722
// helix.toTrackState(alphaCenter, ms.m.p, loc);
683723

684724
TrackState ts_ecal=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, saveTrackStateAtIntercept);
685725
newTrack.getTrackStates().add(ts_ecal);
686-
726+
//if(debug)System.out.println("Helix at ECal from helix.toTrackState");
727+
// if(debug)System.out.println(ts_toTrackState.toString());
728+
// System.out.println("Helix at ECal from helix.toTrackState");
729+
// System.out.println(ts_toTrackState.toString());
730+
731+
// Extrapolate to the ECAL and make a new trackState there.
732+
//BaseTrackState ts_ecal = new BaseTrackState();
733+
//ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber);
734+
687735
Vec targetFace=origin;
688736
Plane targetPlane = new Plane(targetFace, new Vec(0., 1., 0.));
689737
HelixState helixAtTarget=kT.getHelixAtPlane(targetPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane
@@ -695,6 +743,16 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
695743
if(debug)System.out.println("Helix at Target from helix.toTrackState");
696744
if(debug)System.out.println(ts_target.toString());
697745
newTrack.getTrackStates().add(ts_target);
746+
// Extrapolate to Target and make a new trackState there.
747+
/*
748+
BaseTrackState ts_target = new BaseTrackState();
749+
if (target_pos != -999.9 && addTrackStateAtTarget){
750+
ts_target = TrackUtils.getTrackExtrapAtTargetRK(newTrack, target_pos, beamPosition, fM, 0);
751+
if (ts_target != null){
752+
newTrack.getTrackStates().add(ts_target);
753+
}
754+
}
755+
*/
698756

699757
// other track properties
700758
newTrack.setChisq(kT.chi2);

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

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -717,6 +717,14 @@ private void buildH(StateVector S, DMatrixRMaj H) {
717717
}
718718
}
719719

720+
public boolean isInSensor(Double[] rLocal){
721+
double tol = kPar.edgeTolerance; // Tolerance on the check, in mm
722+
// double tol=0.0;
723+
if (rLocal[0] < m.xExtent[0] - tol || rLocal[0] > m.xExtent[1] + tol) return false;
724+
if (rLocal[1] < m.yExtent[0] - tol || rLocal[1] > m.yExtent[1] + tol) return false;
725+
return true;
726+
}
727+
720728
// Comparator functions for sorting measurement sites by layer number
721729
static Comparator<MeasurementSite> SiteComparatorUp = new Comparator<MeasurementSite>() {
722730
public int compare(MeasurementSite s1, MeasurementSite s2) {

0 commit comments

Comments
 (0)