From 093751f4535456f6e909e8ba671fbce54578552b Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 8 May 2025 16:01:12 -0700 Subject: [PATCH 1/4] first set of changes; add switch in propagateRK to return helix at plane intercept; start to use this in TrackStates --- .../hps/recon/tracking/kalman/HelixState.java | 26 +++-- .../hps/recon/tracking/kalman/HelixTest3.java | 2 +- .../hps/recon/tracking/kalman/KalTrack.java | 35 +++++- .../tracking/kalman/KalmanDriverHPS.java | 4 +- .../tracking/kalman/KalmanInterface.java | 103 ++++++++++++++---- .../recon/tracking/kalman/KalmanKinkFit.java | 4 +- .../tracking/kalman/PropagatedTrackState.java | 2 +- 7 files changed, 141 insertions(+), 35 deletions(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java index f26370865c..dead0d178f 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java @@ -322,9 +322,11 @@ static Vec pivotTransform(Vec pivot, Vec a, Vec X0, double alpha, double deltaEo return new Vec(5, aP); } - + HelixState propagateRungeKutta(Plane pln, ArrayList yScat, ArrayList XL, org.lcsim.geometry.FieldMap fM, double [] arcLength){ + return propagateRungeKutta(pln, yScat, XL, fM, arcLength, false); + } // Propagate a helix by Runge-Kutta integration to an arbitrary plane - HelixState propagateRungeKutta(Plane pln, ArrayList yScat, ArrayList XL, org.lcsim.geometry.FieldMap fM, double [] arcLength) { + HelixState propagateRungeKutta(Plane pln, ArrayList yScat, ArrayList XL, org.lcsim.geometry.FieldMap fM, double [] arcLength, boolean pivotAtIntersect) { // pln = plane to where the extrapolation is taking place in global coordinates. // The origin of pln will be the new helix pivot point in global coordinates and the origin of the B-field system. // yScat = input array of y values where scattering in silicon will take place. Only those between the start and finish points @@ -335,7 +337,7 @@ HelixState propagateRungeKutta(Plane pln, ArrayList yScat, ArrayList yScat, ArrayList yScat, ArrayList yScat, ArrayList // In the returned TrackState the reference point gets set to the point on the helix closest to the // original pivot point (e.g. the helix intersection with the plane of silicon). // The pivot of the returned TrackState is always the origin (0,0,0) - TrackState toTrackState(double alphaCenter, Plane pln, int loc) { + TrackState toTrackState(double alphaCenter, Plane pln, int loc, boolean pivotAtIntercept) { + // TrackState toTrackState(double alphaCenter, Plane pln, int loc) { // See TrackState for the different choices for loc (e.g. TrackState.atOther) - return KalmanInterface.toTrackState(this, pln, alphaCenter, loc); + return KalmanInterface.toTrackState(this, pln, alphaCenter, loc, pivotAtIntercept); } // Transform a helix from one pivot to another through a non-uniform B field in several steps diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java index 5b0df88a89..e8989debbc 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java @@ -762,7 +762,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code else if (kF.sites.indexOf(site) == kF.sites.size()-1) { loc = TrackState.AtLastHit; } - TrackState ts = KI.createTrackState(site, loc, true); + TrackState ts = KI.createTrackState(site, loc, true,true); states.add(ts); } TrackState lastState = null; diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java index 74998a3a40..202028d022 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java @@ -20,7 +20,7 @@ import org.hps.util.Pair; import org.lcsim.event.TrackerHit; - +import java.lang.Math; /** * Track followed and fitted by the Kalman filter */ @@ -711,6 +711,39 @@ public double originArcLength() { return arcLength[0]; } + public HelixState getHelixAtPlane(Plane xPlane){ + HelixState helixAtPlane = null; + MeasurementSite closeSite = null; + double delY = 66666.; // y-kalman == z-global + double planeY = xPlane.X().v[1]; + for (MeasurementSite site : SiteList) { + SiModule m = site.m; + // we want closest site to the requested xPlane + if (Math.abs(planeY-m.p.X().v[1])> printExtendedTrackInfo(Track HPStrk) { System.out.format(" Null track state at sensor for this sensor\n"); continue; } - double[] mom = ((BaseTrackState) (tsAtSensor)).computeMomentum(bField); + //MG--5/8/25 computeMomentum doesn't return bfield anymore...use getMomentum() + // double[] mom = ((BaseTrackState) (tsAtSensor)).computeMomentum(bField); + double[] mom = ((BaseTrackState) (tsAtSensor)).getMomentum(); double[] momTransformed = CoordinateTransformations.transformVectorToDetector(new BasicHep3Vector(mom)).v(); Hep3Vector loc = TrackStateUtils.getLocationAtSensor(tsAtSensor, sensor, bField); if (loc == null) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java index 0316e91eaa..805da6f2b5 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java @@ -67,7 +67,7 @@ import org.lcsim.lcio.LCIOConstants; import org.lcsim.recon.cat.util.Const; import org.lcsim.util.Driver; - +import org.hps.recon.tracking.BeamlineConstants; /** * This class provides an interface between hps-java and the Kalman Filter fitting and pattern recognition code. * It can be used to refit the hits on an existing hps track, or it can be used to drive the pattern recognition. @@ -107,7 +107,7 @@ public class KalmanInterface { private static final boolean debug = false; static final double SVTcenter = 505.57; private static final double c = 2.99793e8; // Speed of light in m/s - + static double conFac= 1.0e12 / c; private int runNumber = 14168; public void setRunNumber(int runNumber){ @@ -254,7 +254,6 @@ public KalmanInterface(boolean uniformB, KalmanParams kPar, org.lcsim.geometry.F kPat = new KalmanPatRecHPS(kPar); Vec centerB = KalmanInterface.getField(new Vec(0., SVTcenter, 0.), fM); - double conFac = 1.0e12 / c; alphaCenter = conFac/ centerB.mag(); } @@ -344,7 +343,8 @@ public void clearInterface() { } // Create an HPS TrackState from a Kalman HelixState at the location of a particular SiModule - public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoothed) { + // public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoothed) { + public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoothed, boolean pivotAtIntercept) { // Note that the helix parameters that get stored in the TrackState assume a B-field exactly oriented in the // z direction and a pivot point at the origin (0,0,0). The referencePoint of the TrackState is set to the // intersection point with the detector plane. @@ -357,14 +357,26 @@ public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoot sv = ms.aF; } - return sv.helix.toTrackState(alphaCenter, ms.m.p, loc); + // return sv.helix.toTrackState(alphaCenter, ms.m.p, loc); + return sv.helix.toTrackState(alphaCenter, ms.m.p, loc,pivotAtIntercept); } + public double[][] getCovarianceFromHelix(HelixState helix) { + double[][] M = new double[5][5]; + + for (int i = 0; i < 5; ++i) { + for (int j = 0; j < 5; ++j) { + M[i][j] = helix.C.unsafe_get(i, j); + } + } + return M; + } + // Transform a Kalman helix to an HPS helix rotated to the global frame and with the pivot at the origin // Provide covHPS with 15 elements to get the covariance as well // Provide 3-vector position to get the location in HPS global coordinates of the original helix pivot - static double [] toHPShelix(HelixState helixState, Plane pln, double alphaCenter, double [] covHPS, double[] position) { - final boolean debug = false; + static double [] toHPShelix(HelixState helixState, Plane pln, double alphaCenter, double [] covHPS, double[] position, boolean pivotAtIntercept) { + final boolean debug = true; double phiInt = helixState.planeIntersect(pln); if (Double.isNaN(phiInt)) { Logger logger = Logger.getLogger(KalmanInterface.class.getName()); @@ -431,12 +443,15 @@ public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoot return KalmanInterface.getLCSimParams(finalHelixParams.v, alphaCenter); } - static TrackState toTrackState(HelixState helixState, Plane pln, double alphaCenter, int loc) { + // static TrackState toTrackState(HelixState helixState, Plane pln, double alphaCenter, int loc) { + static TrackState toTrackState(HelixState helixState, Plane pln, double alpha, int loc, boolean pivotAtIntercept) { double [] covHPS = new double[15]; double [] position = new double[3]; - double [] helixHPS = KalmanInterface.toHPShelix(helixState, pln, alphaCenter, covHPS, position); - - return new BaseTrackState( helixHPS, covHPS, position, loc); + double [] helixHPS = KalmanInterface.toHPShelix(helixState, pln, alpha, covHPS, position, pivotAtIntercept); + //position is filled in toHPShelix and used at hte reference in the track state + // return new BaseTrackState( helixHPS, covHPS, position, loc); + double bLocal=conFac/alpha; + return new BaseTrackState( helixHPS, position, covHPS, loc, bLocal); } public void printGBLStripClusterData(GBLStripClusterData clstr) { @@ -575,19 +590,26 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { // and make this trackstate's params the overall track params DMatrixRMaj globalCov = new DMatrixRMaj(kT.originCovariance()); double[] globalParams = kT.originHelixParms(); + Vec origin=kT.helixAtOrigin.origin; // To get the LCSIM curvature parameter, we want the curvature at the center of the SVT (more-or-less the // center of the magnet), so we need to use the magnetic field there to convert from 1/pt to curvature. // Field at the origin => For 2016 this is ~ 0.430612 T // In the center of SVT => For 2016 this is ~ 0.52340 T - Vec Bfield = KalmanInterface.getField(new Vec(0., SVTcenter ,0.), kT.SiteList.get(0).m.Bfield); + // Vec Bfield = KalmanInterface.getField(new Vec(0., SVTcenter ,0.), kT.SiteList.get(0).m.Bfield); + Vec Bfield = KalmanInterface.getField(origin , kT.SiteList.get(0).m.Bfield); double B = Bfield.mag(); - double[] newParams = getLCSimParams(globalParams, alphaCenter); - double[] newCov = getLCSimCov(globalCov, alphaCenter).asPackedArray(true); - TrackState ts = new BaseTrackState(newParams, newCov, new double[]{0., 0., 0.}, TrackState.AtIP); + double alpha = conFac/ B; + // double[] newParams = getLCSimParams(globalParams, alphaCenter); + double[] newParams = getLCSimParams(globalParams, alpha); + // double[] newCov = getLCSimCov(globalCov, alphaCenter).asPackedArray(true); + double[] newCov = getLCSimCov(globalCov, alpha).asPackedArray(true); + double[] refD=new double[]{origin.v[0],origin.v[2],-origin.v[1]}; + TrackState ts = new BaseTrackState(newParams, refD,newCov, TrackState.AtIP); if (ts != null) { newTrack.getTrackStates().add(ts); newTrack.setTrackParameters(ts.getParameters(), B); newTrack.setCovarianceMatrix(new SymmetricMatrix(5, ts.getCovMatrix(), true)); + System.out.println("Track state at perigee : " +ts.toString()); } // Add the hits to the track int firstHit_idx = -1; @@ -640,15 +662,56 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { */ if (loc == TrackState.AtFirstHit || loc == TrackState.AtLastHit || storeTrackStates) { - ts = createTrackState(site, loc, true); + // ts = createTrackState(site, loc, true); + ts = createTrackState(site, loc, true,true); if (ts != null) newTrack.getTrackStates().add(ts); + System.out.println(ts.toString()); } } - + + double zAtEcal = BeamlineConstants.ECAL_FACE; + if (4441 < runNumber && runNumber < 8100) + zAtEcal = BeamlineConstants.ECAL_FACE_ENGINEERING_RUNS; + Vec ecalFace=new Vec(0.,zAtEcal,0.); + Plane ecalPlane = new Plane(ecalFace, new Vec(0., 1., 0.)); + HelixState helixAtEcal=kT.getHelixAtPlane(ecalPlane); + + + // this is how createTrackState (above) goes from measurement to TrackState + // from the measurements. It calls toHPShelix. + // the helix state here must already be propagated + /* helix.toTrackState(alphaCenter, ms.m.p, loc); */ + + // DMatrixRMaj ecalCov = new DMatrixRMaj(kT.ecalCovariance()); + DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal)); + // double[] ecalParams = kT.ecalHelixParms(); + double[] ecalParams = helixAtEcal.a.v.clone(); + Vec BfieldAtEcal = KalmanInterface.getField(helixAtEcal.origin, kT.SiteList.get(0).m.Bfield); + if(debug)System.out.println(this.getClass().getName()+":: helixAtOrigin origin position = "+helixAtEcal.origin.toString()); + double BAtEcal = BfieldAtEcal.mag(); + double alphaAtEcal = conFac/ BAtEcal; + if(debug)System.out.println(this.getClass().getName()+":: B = "+BAtEcal+" alpha= "+alphaAtEcal); + + TrackState ts_toTrackState=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, true); + + + double[] ecalLCSimParams = getLCSimParams(ecalParams, alphaAtEcal); + double[] ecalLCSimCov = getLCSimCov(ecalCov, alphaAtEcal).asPackedArray(true); + double[] refAtEcal = {helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2],helixAtEcal.origin.v[1]}; + TrackState ts_ecal_helix = new BaseTrackState(ecalLCSimParams,refAtEcal, ecalLCSimCov, TrackState.AtCalorimeter, BAtEcal); + if(debug)System.out.println(this.getClass().getName()+":: Uncorrected track state: curvature = "+ts_ecal_helix.getOmega() + +" bField = "+ts_ecal_helix.getBLocal()+" momentum z = "+ts_ecal_helix.getMomentum()[0]); + System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); + System.out.println(ts_ecal_helix.toString()); + newTrack.getTrackStates().add(ts_ecal_helix); + + System.out.println("Helix at ECal from helix.toTrackState"); + System.out.println(ts_toTrackState.toString()); + // Extrapolate to the ECAL and make a new trackState there. - BaseTrackState ts_ecal = new BaseTrackState(); - ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber); - newTrack.getTrackStates().add(ts_ecal); + // BaseTrackState ts_ecal = new BaseTrackState(); + //ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber); + // Extrapolate to Target and make a new trackState there. BaseTrackState ts_target = new BaseTrackState(); diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java index e5c3aea834..4e80561643 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java @@ -242,7 +242,7 @@ public int outerDOF() { return nothing; } MeasurementSite site = innerTrack.sites.get(innerTrack.initialSite); - return KalmanInterface.toHPShelix(site.aS.helix, site.m.p, KI.alphaCenter, null, null); + return KalmanInterface.toHPShelix(site.aS.helix, site.m.p, KI.alphaCenter, null, null,false); } public double [] outerHelix() { // returns a helix with pivot at the origin, in HPS format if (outerTrack == null) { @@ -250,7 +250,7 @@ public int outerDOF() { return nothing; } MeasurementSite site = outerTrack.sites.get(outerTrack.initialSite); - return KalmanInterface.toHPShelix(site.aS.helix, site.m.p, KI.alphaCenter, null, null); + return KalmanInterface.toHPShelix(site.aS.helix, site.m.p, KI.alphaCenter, null, null,false); } public double [] innerMomentum() { if (innerP == null) { diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java index 52d100fd18..ecd38bb4ab 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java @@ -209,7 +209,7 @@ public class PropagatedTrackState { } // Finally turn the HelixState into an HPS TrackState - trackState = newHelixState.toTrackState(alphaCenter, destinationPlane, TrackState.AtOther); + trackState = newHelixState.toTrackState(alphaCenter, destinationPlane, TrackState.AtOther,false); if (debug) printTrackState(trackState,"final"); } From 9768d82d8e3f7999c1a81f16ab01bddf056f5739 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Tue, 13 May 2025 06:49:30 -0700 Subject: [PATCH 2/4] add option to get the helix with pivot point as the plane intersect --- .../hps/recon/tracking/kalman/HelixState.java | 2 +- .../recon/tracking/kalman/KFOutputDriver.java | 218 ++++++++++-------- .../hps/recon/tracking/kalman/KalTrack.java | 10 +- .../tracking/kalman/KalmanInterface.java | 109 +++++---- 4 files changed, 192 insertions(+), 147 deletions(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java index dead0d178f..79f3f71745 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java @@ -337,7 +337,7 @@ HelixState propagateRungeKutta(Plane pln, ArrayList yScat, ArrayList sensorHits) FillKFTrackPlot(trkpFolder+"phi_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.PHI]); //TH2D - Filling - FillKFTrackPlot(trkpFolder+"d0_vs_phi",isTop,charge,trackState.getPhi(),trackState.getD0()); - FillKFTrackPlot(trkpFolder+"d0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getD0()); - FillKFTrackPlot(trkpFolder+"d0_vs_p",isTop,charge,trackp,trackState.getD0()); - - //Ill defined - should be defined wrt bsX and bsY - FillKFTrackPlot(trkpFolder+"d0bs_vs_p",isTop,charge,trackp,helixParametersAtBS[BaseTrack.D0]); - - FillKFTrackPlot(trkpFolder+"z0_vs_p",isTop,charge,trackp,trackState.getZ0()); - FillKFTrackPlot(trkpFolder+"z0bs_vs_p",isTop,charge,trackp,ts_bs.getZ0()); - - //Interesting plot to get a sense where z-vtx is. - //If z0 is referenced to the right BS z location, the slope of vs tanLambda is 0 - FillKFTrackPlot(trkpFolder+"z0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getZ0()); - FillKFTrackPlot(trkpFolder+"z0bs_vs_tanLambda",isTop,charge,trackState.getTanLambda(),ts_bs.getZ0()); - + if(b_doBig2DPlots){ + FillKFTrackPlot(trkpFolder+"d0_vs_phi",isTop,charge,trackState.getPhi(),trackState.getD0()); + FillKFTrackPlot(trkpFolder+"d0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getD0()); + FillKFTrackPlot(trkpFolder+"d0_vs_p",isTop,charge,trackp,trackState.getD0()); + + //Ill defined - should be defined wrt bsX and bsY + FillKFTrackPlot(trkpFolder+"d0bs_vs_p",isTop,charge,trackp,helixParametersAtBS[BaseTrack.D0]); + + FillKFTrackPlot(trkpFolder+"z0_vs_p",isTop,charge,trackp,trackState.getZ0()); + FillKFTrackPlot(trkpFolder+"z0bs_vs_p",isTop,charge,trackp,ts_bs.getZ0()); + + //Interesting plot to get a sense where z-vtx is. + //If z0 is referenced to the right BS z location, the slope of vs tanLambda is 0 + FillKFTrackPlot(trkpFolder+"z0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getZ0()); + FillKFTrackPlot(trkpFolder+"z0bs_vs_tanLambda",isTop,charge,trackState.getTanLambda(),ts_bs.getZ0()); + } if(b_doRawHitPlots){ @@ -746,8 +748,10 @@ private void doKFresiduals(Track trk, Map sensorHits, E // after the transformation x and y in the sensor frame are reversed // This plot is ill defined. - aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitTrackPosSensor.y(), hitTrackPosSensor.x()); - //aidaKF.histogram2D("hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitPos.y(), hitPos.x()); + if(b_doBig2DPlots){ + aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitTrackPosSensor.y(), hitTrackPosSensor.x()); + } + //aidaKF.histogram2D("hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitPos.y(), hitPos.x()); //aidaKF.histogram2D("hit y vs x lab-frame " + sensor.getName()).fill(hitPos.y(), hitPos.x()); @@ -760,23 +764,28 @@ private void doKFresiduals(Track trk, Map sensorHits, E extrapPosSensor = new BasicHep3Vector(extrapPos.v()); trans.transform(extrapPosSensor); //aidaKF.histogram2D("residual after KF vs u predicted " + sensor.getName()).fill(extrapPosSensor.x(), res); - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); - // select track charge - if(trk.getCharge()>0) { - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); - }else if(trk.getCharge()<0) { - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); - } - + if(b_doBig2DPlots){ + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + } + // select track charge + if(b_doBig2DPlots){ + if(trk.getCharge()>0) { + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + }else if(trk.getCharge()<0) { + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + } + } // post-KF residual Hep3Vector hitPos = new BasicHep3Vector(sensorHits.get(sensor).getPosition()); Hep3Vector hitPosSensor = new BasicHep3Vector(hitPos.v()); trans.transform(hitPosSensor); Hep3Vector resSensor = VecOp.sub(hitPosSensor, extrapPosSensor); - aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName()).fill(extrapPosSensor.y(), resSensor.x()); - aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName()).fill(hitPosSensor.x(), resSensor.x()); - aidaKF.histogram1D(resFolder+"residual_after_KF_" + sensor.getName()).fill(resSensor.x()); - + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName()).fill(extrapPosSensor.y(), resSensor.x()); + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName()).fill(hitPosSensor.x(), resSensor.x()); + } + aidaKF.histogram1D(resFolder+"residual_after_KF_" + sensor.getName()).fill(resSensor.x()); + if (debug) { @@ -888,12 +897,13 @@ private void doKFresiduals(Track trk, Map sensorHits, E int spacing = 0; if (vol == "_bottom") spacing = sensors.size()/2 + mod; - - aidaKF.histogram2D(resFolder+"uresidual_KF_mod").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); - aidaKF.profile1D(resFolder+"uresidual_KF_mod_p").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); - aidaKF.histogram1D(resFolder+"uresidual_KF_" + sensorName).fill(trackRes.getDoubleVal(i_hit)); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensorName).fill(hitPosSensorG.x(),trackRes.getDoubleVal(i_hit)); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensorName).fill(extrapPosSensor.y(),trackRes.getDoubleVal(i_hit)); + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"uresidual_KF_mod").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensorName).fill(hitPosSensorG.x(),trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensorName).fill(extrapPosSensor.y(),trackRes.getDoubleVal(i_hit)); + } + aidaKF.histogram1D(resFolder+"uresidual_KF_" + sensorName).fill(trackRes.getDoubleVal(i_hit)); + aidaKF.profile1D(resFolder+"uresidual_KF_mod_p").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); aidaKF.histogram1D(epullFolder+"ureserror_KF_" + sensorName).fill(trackRes.getFloatVal(i_hit)); aidaKF.histogram1D(epullFolder+"ures_pull_KF_" + sensorName).fill(trackRes.getDoubleVal(i_hit) / Math.sqrt(trackRes.getFloatVal(i_hit))); @@ -904,10 +914,10 @@ private void doKFresiduals(Track trk, Map sensorHits, E double dT_hit_track = hitTime - trackTime; double dT_hit_sigma = (hitTime - trackTime) / trackTimeSD; - - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_"+sensorName).fill(dT_hit_track,trackRes.getDoubleVal(i_hit)); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_"+sensorName).fill(dT_hit_sigma,trackRes.getDoubleVal(i_hit)); - + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_"+sensorName).fill(dT_hit_track,trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_"+sensorName).fill(dT_hit_sigma,trackRes.getDoubleVal(i_hit)); + } @@ -973,7 +983,7 @@ private void setupEoPPlots() { aidaKF.histogram2D(eopFolder+"EoP_vs_trackP"+charge+vol,200,0,6,200,0,2); aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda"+charge+vol,200,lmin,lmax,200,0,2); aidaKF.histogram2D(eopFolder+"EoP_vs_phi"+charge+vol,200,-0.2,0.2,200,0,2); - } + } aidaKF.histogram1D(eopFolder+"Ecluster"+vol+"_fid",200,0,5); aidaKF.histogram1D(eopFolder+"EoP"+vol+"_fid",200,0,2); @@ -1072,7 +1082,9 @@ private void setupPlots() { // aidaKF.histogram2D(resFolder+"bresidual_KF_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5, nbins, -xmax,xmax); // aidaKF.profile1D(resFolder+"bresidual_KF_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); - aidaKF.histogram2D(resFolder+"uresidual_KF_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5, 400, -0.4,0.4); + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"uresidual_KF_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5, 400, -0.4,0.4); + } aidaKF.profile1D(resFolder+"uresidual_KF_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); @@ -1113,24 +1125,25 @@ private void setupPlots() { aidaKF.histogram1D(resFolder+"residual_after_KF_" + sensor.getName(), nbins, -xmax, xmax); // aidaKF.histogram1D(resFolder+"bresidual_KF_" + sensor.getName(), nbins, -xmax, xmax); aidaKF.histogram1D(resFolder+"uresidual_KF_" + sensor.getName(), nbins, -xmax, xmax); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensor.getName(),100,-20.0,20.0,100,-0.1,0.1); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensor.getName(),300,-60.0,60.0,100,-0.1,0.1); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_" + sensor.getName(),100,-10.0,10.0,100,-0.1,0.1); - aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_" + sensor.getName(),100,-5.0,5.0,100,-0.1,0.1); - + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensor.getName(),100,-20.0,20.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensor.getName(),300,-60.0,60.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_" + sensor.getName(),100,-10.0,10.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_" + sensor.getName(),100,-5.0,5.0,100,-0.1,0.1); + } // aidaKF.histogram1D(epullFolder+"breserror_KF_" + sensor.getName(), nbins, 0.0, 0.1); aidaKF.histogram1D(epullFolder+"ureserror_KF_" + sensor.getName(), nbins, 0.0, 0.2); // aidaKF.histogram1D(epullFolder+"bres_pull_KF_" + sensor.getName(), nbins, -5, 5); aidaKF.histogram1D(epullFolder+"ures_pull_KF_" + sensor.getName(), nbins, -5, 5); - - aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName(), 100, -20.0, 20.0, 100, -0.04, 0.04); - aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName(), 100, -55.0, 55.0, 100, -0.04, 0.04); - aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName(), 300, -60.0, 60.0, 300, -25, 25); - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); - aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); - + if(b_doBig2DPlots){ + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName(), 100, -20.0, 20.0, 100, -0.04, 0.04); + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName(), 100, -55.0, 55.0, 100, -0.04, 0.04); + aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName(), 300, -60.0, 60.0, 300, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + } aidaKF.histogram1D(hitFolder+"raw_hit_t0_"+sensor.getName(),200, -100, 100.0); aidaKF.histogram1D(hitFolder+"raw_hit_amplitude_"+sensor.getName(),200, 0.0, 4000.0); aidaKF.histogram1D(hitFolder+"raw_hit_chisq_"+sensor.getName(),200, 0.0, 2.0); @@ -1214,25 +1227,26 @@ private void setupPlots() { //TH2Ds - aidaKF.histogram2D(trkpFolder+"d0_vs_phi"+vol+charge,nbins_t,-0.3,0.3,nbins_t,-5.0,5.0); - aidaKF.histogram2D(trkpFolder+"Chi2_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t*2,0,200); - //aidaKF.histogram2D("d0_vs_phi_bs"+vol+charge,nbins_t,-5.0,5.0,nbins_t,-0.3,0.3); - aidaKF.histogram2D(trkpFolder+"d0_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_t,-5.0,5.0); - aidaKF.histogram2D(trkpFolder+"d0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); - aidaKF.histogram2D(trkpFolder+"d0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-5.0,5.0); - aidaKF.histogram2D(trkpFolder+"z0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); - aidaKF.histogram2D(trkpFolder+"z0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-z0bsmax,z0bsmax); - aidaKF.histogram2D(trkpFolder+"z0_vs_tanLambda"+vol+charge, nbins_t,-0.1,0.1,nbins_t,-z0max,z0max); - aidaKF.histogram2D(trkpFolder+"z0bs_vs_tanLambda"+vol+charge,nbins_t,-0.1,0.1,nbins_t,-z0bsmax,z0bsmax); - - aidaKF.histogram2D(trkpFolder+"p_Missing1Hit"+vol+charge,8,0,8,nbins_p,0.0,pmax); - aidaKF.histogram2D(trkpFolder+"p_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); - aidaKF.histogram2D(trkpFolder+"p_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); - // aidaKF.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); - - aidaKF.histogram2D(trkpFolder+"pT_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); - aidaKF.histogram2D(trkpFolder+"pT_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); - + if(b_doBig2DPlots){ + aidaKF.histogram2D(trkpFolder+"d0_vs_phi"+vol+charge,nbins_t,-0.3,0.3,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"Chi2_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t*2,0,200); + //aidaKF.histogram2D("d0_vs_phi_bs"+vol+charge,nbins_t,-5.0,5.0,nbins_t,-0.3,0.3); + aidaKF.histogram2D(trkpFolder+"d0_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"d0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"d0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"z0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"z0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-z0bsmax,z0bsmax); + aidaKF.histogram2D(trkpFolder+"z0_vs_tanLambda"+vol+charge, nbins_t,-0.1,0.1,nbins_t,-z0max,z0max); + aidaKF.histogram2D(trkpFolder+"z0bs_vs_tanLambda"+vol+charge,nbins_t,-0.1,0.1,nbins_t,-z0bsmax,z0bsmax); + + aidaKF.histogram2D(trkpFolder+"p_Missing1Hit"+vol+charge,8,0,8,nbins_p,0.0,pmax); + aidaKF.histogram2D(trkpFolder+"p_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaKF.histogram2D(trkpFolder+"p_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + // aidaKF.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); + + aidaKF.histogram2D(trkpFolder+"pT_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaKF.histogram2D(trkpFolder+"pT_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + } if (b_doDetailPlots) { diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java index 202028d022..72baaf0b80 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java @@ -710,8 +710,8 @@ public double originArcLength() { } return arcLength[0]; } - - public HelixState getHelixAtPlane(Plane xPlane){ + + public HelixState getHelixAtPlane(Plane xPlane, boolean pivotAtIntersect){ HelixState helixAtPlane = null; MeasurementSite closeSite = null; double delY = 66666.; // y-kalman == z-global @@ -734,7 +734,7 @@ public HelixState getHelixAtPlane(Plane xPlane){ } } double[] arcL = new double[1]; - helixAtPlane = closeSite.aS.helix.propagateRungeKutta(xPlane, yScat, XLscat, closeSite.m.Bfield, arcL); + helixAtPlane = closeSite.aS.helix.propagateRungeKutta(xPlane, yScat, XLscat, closeSite.m.Bfield, arcL, pivotAtIntersect); if (debug) { System.out.format("KalTrack::getHelixAtPlane: arc length to the first measurement = %9.4f\n", arcL[0]); } @@ -1673,4 +1673,8 @@ public boolean equals(Object other) { // Consider two tracks to be equal if they public int hashCode() { return nHits + 100 * ID; } + + public KalmanParams getKalPars(){ + return kPar; + } } diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java index 805da6f2b5..830f2597a3 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java @@ -109,6 +109,7 @@ public class KalmanInterface { private static final double c = 2.99793e8; // Speed of light in m/s static double conFac= 1.0e12 / c; private int runNumber = 14168; + private boolean saveTrackStateAtIntercept = true; public void setRunNumber(int runNumber){ this.runNumber = runNumber; @@ -310,11 +311,22 @@ public static Vec localHpsToKal(double[] HPSvec) { public static double[] localKalToHps(Vec KalVec) { double[] HPSvec = new double[3]; HPSvec[0] = KalVec.v[1]; - HPSvec[1] = -KalVec.v[0]; - HPSvec[2] = KalVec.v[2]; + // HPSvec[1] = -KalVec.v[0]; + //HPSvec[2] = KalVec.v[2]; + // think the sign above is switched + HPSvec[1] = KalVec.v[0]; + HPSvec[2] = -KalVec.v[2]; return HPSvec; + } + + public static double[] hpsGlobalToTracking(double[] hpsGbl){ + double[] hpsTrk = new double[3]; + hpsTrk[0] = hpsGbl[2]; + hpsTrk[1] = hpsGbl[0]; + hpsTrk[2] = hpsGbl[1]; + return hpsTrk; } - + // Return the entire list of Kalman SiModule public ArrayList getSiModuleList() { return SiMlist; @@ -356,9 +368,9 @@ public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoot if (!ms.filtered) return null; sv = ms.aF; } - - // return sv.helix.toTrackState(alphaCenter, ms.m.p, loc); - return sv.helix.toTrackState(alphaCenter, ms.m.p, loc,pivotAtIntercept); + Vec planeCenter=ms.m.p.X(); + double alpha=conFac/KalmanInterface.getField(ms.m.p.X(),ms.m.Bfield).mag(); + return sv.helix.toTrackState(alpha, ms.m.p, loc,pivotAtIntercept); } public double[][] getCovarianceFromHelix(HelixState helix) { @@ -375,8 +387,10 @@ public double[][] getCovarianceFromHelix(HelixState helix) { // Transform a Kalman helix to an HPS helix rotated to the global frame and with the pivot at the origin // Provide covHPS with 15 elements to get the covariance as well // Provide 3-vector position to get the location in HPS global coordinates of the original helix pivot - static double [] toHPShelix(HelixState helixState, Plane pln, double alphaCenter, double [] covHPS, double[] position, boolean pivotAtIntercept) { - final boolean debug = true; + // 5/10/25: add option to return hps helix with reference at plane intercept -- mgraham + // 5/10/25: change variable name from alphaCenter to alpha to avoid confusion -- mgraham + static double [] toHPShelix(HelixState helixState, Plane pln, double alpha, double [] covHPS, double[] position, boolean pivotAtIntercept) { + final boolean debug = false; double phiInt = helixState.planeIntersect(pln); if (Double.isNaN(phiInt)) { Logger logger = Logger.getLogger(KalmanInterface.class.getName()); @@ -417,30 +431,35 @@ public double[][] getCovarianceFromHelix(HelixState helix) { // Pivot transform to the final pivot at the origin Vec finalPivot = new Vec(0.,0.,0.); - Vec finalHelixParams = HelixState.pivotTransform(finalPivot, helixParamsRotated, pivotGlobal, alphaCenter, 0.); - HelixState.makeF(finalHelixParams, F, helixParamsRotated, alphaCenter); + Vec finalHelixParams = null; + if(pivotAtIntercept) + finalHelixParams = helixParamsRotated; + else + finalHelixParams = HelixState.pivotTransform(finalPivot, helixParamsRotated, pivotGlobal, alpha, 0.); + + HelixState.makeF(finalHelixParams, F, helixParamsRotated, alpha); CommonOps_DDRM.multTransB(covRotated, F, tempM); CommonOps_DDRM.mult(F, tempM, covRotated); if (debug) { finalPivot.print("final pivot point"); finalHelixParams.print("final helix parameters"); HelixPlaneIntersect hpi = new HelixPlaneIntersect(); - phiInt = hpi.planeIntersect(finalHelixParams, finalPivot, alphaCenter, pln); + phiInt = hpi.planeIntersect(finalHelixParams, finalPivot, alpha, pln); if (!Double.isNaN(phiInt)) { - Vec rInt = HelixState.atPhi(finalPivot, finalHelixParams, phiInt, alphaCenter); + Vec rInt = HelixState.atPhi(finalPivot, finalHelixParams, phiInt, alpha); rInt.print("final helix intersection with given plane"); } System.out.format("Exiting KalmanInterface.toHPShelix\n"); } if (covHPS != null) { - double [] temp = KalmanInterface.getLCSimCov(covRotated, alphaCenter).asPackedArray(true); + double [] temp = KalmanInterface.getLCSimCov(covRotated, alpha).asPackedArray(true); for (int i=0; i<15; ++i) covHPS[i] = temp[i]; } if (position != null) { double [] temp = KalmanInterface.vectorKalmanToGlb(pivotGlobal); for (int i=0; i<3; ++i) position[i] = temp[i]; } - return KalmanInterface.getLCSimParams(finalHelixParams.v, alphaCenter); + return KalmanInterface.getLCSimParams(finalHelixParams.v, alpha); } // static TrackState toTrackState(HelixState helixState, Plane pln, double alphaCenter, int loc) { @@ -448,10 +467,11 @@ static TrackState toTrackState(HelixState helixState, Plane pln, double alpha, i double [] covHPS = new double[15]; double [] position = new double[3]; double [] helixHPS = KalmanInterface.toHPShelix(helixState, pln, alpha, covHPS, position, pivotAtIntercept); - //position is filled in toHPShelix and used at hte reference in the track state - // return new BaseTrackState( helixHPS, covHPS, position, loc); + //position is filled in toHPShelix and used at the reference in the track state + //however it is in the global frame and is expected in the HPS tracking, so fix that + double[] posTrking=hpsGlobalToTracking(position); double bLocal=conFac/alpha; - return new BaseTrackState( helixHPS, position, covHPS, loc, bLocal); + return new BaseTrackState( helixHPS, posTrking, covHPS, loc, bLocal); } public void printGBLStripClusterData(GBLStripClusterData clstr) { @@ -596,20 +616,23 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { // Field at the origin => For 2016 this is ~ 0.430612 T // In the center of SVT => For 2016 this is ~ 0.52340 T // Vec Bfield = KalmanInterface.getField(new Vec(0., SVTcenter ,0.), kT.SiteList.get(0).m.Bfield); - Vec Bfield = KalmanInterface.getField(origin , kT.SiteList.get(0).m.Bfield); + // double[] refD=new double[]{origin.v[0],-origin.v[2],origin.v[1]}; + double[] refD=localKalToHps(origin); + Vec originHPS=new Vec(3,refD); + Vec Bfield = KalmanInterface.getField(originHPS, kT.SiteList.get(0).m.Bfield); double B = Bfield.mag(); double alpha = conFac/ B; // double[] newParams = getLCSimParams(globalParams, alphaCenter); double[] newParams = getLCSimParams(globalParams, alpha); // double[] newCov = getLCSimCov(globalCov, alphaCenter).asPackedArray(true); double[] newCov = getLCSimCov(globalCov, alpha).asPackedArray(true); - double[] refD=new double[]{origin.v[0],origin.v[2],-origin.v[1]}; - TrackState ts = new BaseTrackState(newParams, refD,newCov, TrackState.AtIP); + // TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtIP); + TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtIP, B); if (ts != null) { newTrack.getTrackStates().add(ts); newTrack.setTrackParameters(ts.getParameters(), B); newTrack.setCovarianceMatrix(new SymmetricMatrix(5, ts.getCovMatrix(), true)); - System.out.println("Track state at perigee : " +ts.toString()); + if(debug)System.out.println("Track state at perigee : " +ts.toString()); } // Add the hits to the track int firstHit_idx = -1; @@ -662,10 +685,9 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { */ if (loc == TrackState.AtFirstHit || loc == TrackState.AtLastHit || storeTrackStates) { - // ts = createTrackState(site, loc, true); - ts = createTrackState(site, loc, true,true); + ts = createTrackState(site, loc, true, saveTrackStateAtIntercept); if (ts != null) newTrack.getTrackStates().add(ts); - System.out.println(ts.toString()); + if(debug)System.out.println(ts.toString()); } } @@ -674,39 +696,43 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { zAtEcal = BeamlineConstants.ECAL_FACE_ENGINEERING_RUNS; Vec ecalFace=new Vec(0.,zAtEcal,0.); Plane ecalPlane = new Plane(ecalFace, new Vec(0., 1., 0.)); - HelixState helixAtEcal=kT.getHelixAtPlane(ecalPlane); - - - // this is how createTrackState (above) goes from measurement to TrackState - // from the measurements. It calls toHPShelix. - // the helix state here must already be propagated - /* helix.toTrackState(alphaCenter, ms.m.p, loc); */ + HelixState helixAtEcal=kT.getHelixAtPlane(ecalPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane // DMatrixRMaj ecalCov = new DMatrixRMaj(kT.ecalCovariance()); DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal)); // double[] ecalParams = kT.ecalHelixParms(); double[] ecalParams = helixAtEcal.a.v.clone(); Vec BfieldAtEcal = KalmanInterface.getField(helixAtEcal.origin, kT.SiteList.get(0).m.Bfield); - if(debug)System.out.println(this.getClass().getName()+":: helixAtOrigin origin position = "+helixAtEcal.origin.toString()); + if(debug)System.out.println(this.getClass().getName()+":: helixAtOrigin origin position @ ecal = "+helixAtEcal.origin.toString()); double BAtEcal = BfieldAtEcal.mag(); double alphaAtEcal = conFac/ BAtEcal; - if(debug)System.out.println(this.getClass().getName()+":: B = "+BAtEcal+" alpha= "+alphaAtEcal); - - TrackState ts_toTrackState=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, true); - + if(debug)System.out.println(this.getClass().getName()+":: B = "+BAtEcal+" alpha= "+alphaAtEcal); double[] ecalLCSimParams = getLCSimParams(ecalParams, alphaAtEcal); double[] ecalLCSimCov = getLCSimCov(ecalCov, alphaAtEcal).asPackedArray(true); - double[] refAtEcal = {helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2],helixAtEcal.origin.v[1]}; + // double[] refAtEcal = {helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2],helixAtEcal.origin.v[1]}; + // double[] refAtEcal = {helixAtEcal.origin.v[1],helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2]}; + double[] refAtEcal = localKalToHps(helixAtEcal.origin); + if(debug)System.out.println(this.getClass().getName()+":: reference to TrackState @ ecal = "+ + +refAtEcal[0]+", "+refAtEcal[1]+","+ refAtEcal[2]); TrackState ts_ecal_helix = new BaseTrackState(ecalLCSimParams,refAtEcal, ecalLCSimCov, TrackState.AtCalorimeter, BAtEcal); if(debug)System.out.println(this.getClass().getName()+":: Uncorrected track state: curvature = "+ts_ecal_helix.getOmega() +" bField = "+ts_ecal_helix.getBLocal()+" momentum z = "+ts_ecal_helix.getMomentum()[0]); - System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); - System.out.println(ts_ecal_helix.toString()); + if(debug)System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); + if(debug)System.out.println(ts_ecal_helix.toString()); + if(debug)System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); + if(debug)System.out.println(ts_ecal_helix.toString()); newTrack.getTrackStates().add(ts_ecal_helix); - System.out.println("Helix at ECal from helix.toTrackState"); - System.out.println(ts_toTrackState.toString()); + + // this is how createTrackState (above) goes from measurement to TrackState + // from the measurements. It calls toHPShelix. + // the helix state here must already be propagated + /* helix.toTrackState(alphaCenter, ms.m.p, loc); */ + + //TrackState ts_toTrackState=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, saveTrackStateAtIntercept); + //if(debug)System.out.println("Helix at ECal from helix.toTrackState"); + //if(debug)System.out.println(ts_toTrackState.toString()); // Extrapolate to the ECAL and make a new trackState there. // BaseTrackState ts_ecal = new BaseTrackState(); @@ -822,6 +848,7 @@ private void addHitsToTrack(BaseTrack newTrack) { } // Create an HPS track from a Kalman seed + // mg 5/12/2025 ... I'm not sure this is used anywhere so keep alphaCenter for now public BaseTrack createTrack(SeedTrack trk) { double[] newPivot = { 0., 0., 0. }; double[] params = getLCSimParams(trk.pivotTransform(newPivot), alphaCenter); From 9cccbc330bdf7f006b51b4481c8425a32ecb0533 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 21 May 2025 17:23:22 -0700 Subject: [PATCH 3/4] add blocal and make computeMomentum use only class variables --- .../examples/StarterAnalysisDriver.java | 2 +- .../recon/tracking/kalman/KFOutputDriver.java | 20 ++++-- .../tracking/kalman/KalmanInterface.java | 68 ++++++++++--------- 3 files changed, 50 insertions(+), 40 deletions(-) diff --git a/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java b/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java index 0018c6a072..337593eb9c 100644 --- a/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java +++ b/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java @@ -147,7 +147,7 @@ private List processTracks(EventHeader event) { for (Track track : tracklist) { //remember, these tracks are in the lcsim tracking frame! BaseTrackState ts = (BaseTrackState) track.getTrackStates().get(0); - ts.computeMomentum(bfield); + ts.computeMomentum(); tkanalList.add(new LCIOTrackAnalysis(track, hittomc, hittostrip, hittorotated)); } return tkanalList; diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java index ee6391b81e..a3ec010ecf 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java @@ -78,9 +78,9 @@ public class KFOutputDriver extends Driver { String hitFolder="/hit/"; String eopFolder = "/EoP/"; // private boolean b_doKFkinks = false; - private boolean b_doKFresiduals = true; + private boolean b_doKFresiduals = false; private boolean b_doDetailPlots = false; - private boolean b_doRawHitPlots = true; + private boolean b_doRawHitPlots = false; private boolean b_doBig2DPlots = false; //The field map for extrapolation private FieldMap bFieldMap; @@ -270,7 +270,9 @@ public void process(EventHeader event) { if (particle.getTracks().isEmpty() || particle.getClusters().isEmpty()) continue; Track track = particle.getTracks().get(0); - Cluster cluster = particle.getClusters().get(0); + Cluster cluster = particle.getClusters().get(0); + if(debug) + System.out.println(this.getClass().getName()+":: adding track and cluster to lists"); tracks.add(track); TrackClusterPairs.put(track,cluster); } @@ -284,7 +286,8 @@ public void process(EventHeader event) { RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event); for (Track trk : tracks) { - + if(debug) System.out.println(this.getClass().getName()+":: track chi2 = "+trk.getChi2()); + if(debug) System.out.println(this.getClass().getName()+":: track nHits = "+trk.getTrackerHits().size()); if (trk.getChi2() > chi2Cut) continue; @@ -294,8 +297,11 @@ public void process(EventHeader event) { if(debug) System.out.println("Track passed hits d0 = "+trk.getTrackStates().get(0).getD0()); - - Hep3Vector momentum = new BasicHep3Vector(trk.getTrackStates().get(0).getMomentum()); + + if(debug)System.out.println(this.getClass().getName()+":: local B field = "+trk.getTrackStates().get(0).getBLocal()); + // ((BaseTrackState)trk.getTrackStates().get(0)).computeMomentum(trk.getTrackStates().get(0).getBLocal()); + Hep3Vector momentum = new BasicHep3Vector(trk.getTrackStates().get(0).getMomentum()); + if(debug) System.out.println(this.getClass().getName()+":: track momentum = "+momentum.magnitude()); if (momentum.magnitude() < minMom) continue; @@ -337,7 +343,7 @@ public void process(EventHeader event) { System.out.printf("TrackerHit null sensor %s \n", hit.toString()); } _trkTimeSigma=getTrackTime(sensorHits); - doBasicKFtrack(trk,sensorHits); + doBasicKFtrack(trk,sensorHits); if (b_doKFresiduals) doKFresiduals(trk, sensorHits,event); diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java index 830f2597a3..02acbedcd0 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java @@ -606,28 +606,19 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { kT.sortSites(true); BaseTrack newTrack = new BaseTrack(); - // Add trackstate at IP as first trackstate, + // Add perigee trackstate at IP as first trackstate, // and make this trackstate's params the overall track params DMatrixRMaj globalCov = new DMatrixRMaj(kT.originCovariance()); double[] globalParams = kT.originHelixParms(); - Vec origin=kT.helixAtOrigin.origin; - // To get the LCSIM curvature parameter, we want the curvature at the center of the SVT (more-or-less the - // center of the magnet), so we need to use the magnetic field there to convert from 1/pt to curvature. - // Field at the origin => For 2016 this is ~ 0.430612 T - // In the center of SVT => For 2016 this is ~ 0.52340 T - // Vec Bfield = KalmanInterface.getField(new Vec(0., SVTcenter ,0.), kT.SiteList.get(0).m.Bfield); - // double[] refD=new double[]{origin.v[0],-origin.v[2],origin.v[1]}; + Vec origin=kT.helixAtOrigin.origin; double[] refD=localKalToHps(origin); Vec originHPS=new Vec(3,refD); Vec Bfield = KalmanInterface.getField(originHPS, kT.SiteList.get(0).m.Bfield); double B = Bfield.mag(); double alpha = conFac/ B; - // double[] newParams = getLCSimParams(globalParams, alphaCenter); double[] newParams = getLCSimParams(globalParams, alpha); - // double[] newCov = getLCSimCov(globalCov, alphaCenter).asPackedArray(true); double[] newCov = getLCSimCov(globalCov, alpha).asPackedArray(true); - // TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtIP); - TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtIP, B); + TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtPerigee, B); if (ts != null) { newTrack.getTrackStates().add(ts); newTrack.setTrackParameters(ts.getParameters(), B); @@ -691,55 +682,67 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { } } + /* define ecal plane, getHelix, and b-field conversion value */ double zAtEcal = BeamlineConstants.ECAL_FACE; if (4441 < runNumber && runNumber < 8100) zAtEcal = BeamlineConstants.ECAL_FACE_ENGINEERING_RUNS; Vec ecalFace=new Vec(0.,zAtEcal,0.); Plane ecalPlane = new Plane(ecalFace, new Vec(0., 1., 0.)); HelixState helixAtEcal=kT.getHelixAtPlane(ecalPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane - - // DMatrixRMaj ecalCov = new DMatrixRMaj(kT.ecalCovariance()); - DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal)); - // double[] ecalParams = kT.ecalHelixParms(); - double[] ecalParams = helixAtEcal.a.v.clone(); Vec BfieldAtEcal = KalmanInterface.getField(helixAtEcal.origin, kT.SiteList.get(0).m.Bfield); if(debug)System.out.println(this.getClass().getName()+":: helixAtOrigin origin position @ ecal = "+helixAtEcal.origin.toString()); double BAtEcal = BfieldAtEcal.mag(); double alphaAtEcal = conFac/ BAtEcal; - if(debug)System.out.println(this.getClass().getName()+":: B = "+BAtEcal+" alpha= "+alphaAtEcal); + + /* + DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal)); + double[] ecalParams = helixAtEcal.a.v.clone(); double[] ecalLCSimParams = getLCSimParams(ecalParams, alphaAtEcal); double[] ecalLCSimCov = getLCSimCov(ecalCov, alphaAtEcal).asPackedArray(true); - // double[] refAtEcal = {helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2],helixAtEcal.origin.v[1]}; - // double[] refAtEcal = {helixAtEcal.origin.v[1],helixAtEcal.origin.v[0],-helixAtEcal.origin.v[2]}; double[] refAtEcal = localKalToHps(helixAtEcal.origin); if(debug)System.out.println(this.getClass().getName()+":: reference to TrackState @ ecal = "+ +refAtEcal[0]+", "+refAtEcal[1]+","+ refAtEcal[2]); TrackState ts_ecal_helix = new BaseTrackState(ecalLCSimParams,refAtEcal, ecalLCSimCov, TrackState.AtCalorimeter, BAtEcal); if(debug)System.out.println(this.getClass().getName()+":: Uncorrected track state: curvature = "+ts_ecal_helix.getOmega() +" bField = "+ts_ecal_helix.getBLocal()+" momentum z = "+ts_ecal_helix.getMomentum()[0]); - if(debug)System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); - if(debug)System.out.println(ts_ecal_helix.toString()); - if(debug)System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); - if(debug)System.out.println(ts_ecal_helix.toString()); - newTrack.getTrackStates().add(ts_ecal_helix); - + System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); + System.out.println(ts_ecal_helix.toString()); + */ + //newTrack.getTrackStates().add(ts_ecal_helix); + // this is how createTrackState (above) goes from measurement to TrackState // from the measurements. It calls toHPShelix. // the helix state here must already be propagated - /* helix.toTrackState(alphaCenter, ms.m.p, loc); */ + // helix.toTrackState(alphaCenter, ms.m.p, loc); - //TrackState ts_toTrackState=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, saveTrackStateAtIntercept); + TrackState ts_ecal=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, saveTrackStateAtIntercept); + newTrack.getTrackStates().add(ts_ecal); //if(debug)System.out.println("Helix at ECal from helix.toTrackState"); - //if(debug)System.out.println(ts_toTrackState.toString()); + // if(debug)System.out.println(ts_toTrackState.toString()); + // System.out.println("Helix at ECal from helix.toTrackState"); + // System.out.println(ts_toTrackState.toString()); // Extrapolate to the ECAL and make a new trackState there. - // BaseTrackState ts_ecal = new BaseTrackState(); + //BaseTrackState ts_ecal = new BaseTrackState(); //ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber); - - + + Vec targetFace=origin; + Plane targetPlane = new Plane(targetFace, new Vec(0., 1., 0.)); + targetFace.print("target face"); + HelixState helixAtTarget=kT.getHelixAtPlane(targetPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane + helixAtTarget.origin.print("helixAtTarget.origin"); + Vec BfieldAtTarget = KalmanInterface.getField(helixAtTarget.origin, kT.SiteList.get(0).m.Bfield); + if(debug)System.out.println(this.getClass().getName()+":: helixAtOrigin origin position @ target = "+helixAtTarget.origin.toString()); + double BAtTarget = BfieldAtTarget.mag(); + double alphaAtTarget = conFac/ BAtTarget; + TrackState ts_target=helixAtTarget.toTrackState(alphaAtTarget, targetPlane, TrackState.AtTarget, saveTrackStateAtIntercept); + System.out.println("Helix at Target from helix.toTrackState"); + System.out.println(ts_target.toString()); + newTrack.getTrackStates().add(ts_target); // Extrapolate to Target and make a new trackState there. + /* BaseTrackState ts_target = new BaseTrackState(); if (target_pos != -999.9 && addTrackStateAtTarget){ ts_target = TrackUtils.getTrackExtrapAtTargetRK(newTrack, target_pos, beamPosition, fM, 0); @@ -747,6 +750,7 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { newTrack.getTrackStates().add(ts_target); } } + */ // other track properties newTrack.setChisq(kT.chi2); From 62a880aab57c37f63fda4faed2eca6987858b31f Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 30 May 2025 09:49:40 -0700 Subject: [PATCH 4/4] use bLocal for the vertexing --- .../hps/recon/particle/HpsReconParticleDriver.java | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java b/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java index 54c93b14a6..b562e5896a 100644 --- a/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java +++ b/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java @@ -509,12 +509,14 @@ private BilliorVertex fitVertex(Constraint constraint, ReconstructedParticle ele BilliorTrack electronBTrack = toBilliorTrack(electron.getTracks().get(0)); BilliorTrack positronBTrack = toBilliorTrack(positron.getTracks().get(0)); - // Create a vertex fitter from the magnetic field. + // Create a vertex fitter from the magnetic field. // Note that the vertexing code uses the tracking frame coordinates // HPS X => TRACK Y // HPS Y => TRACK Z // HPS Z => TRACK X - BilliorVertexer vtxFitter = new BilliorVertexer(bField); + //first get the field @ perigee reference from the first tracks state (doesn't matter which one) + double bLocal=(electron.getTracks().get(0)).getTrackStates().get(TrackState.AtPerigee).getBLocal(); + BilliorVertexer vtxFitter = new BilliorVertexer(bLocal); // TODO: The beam size should come from the conditions database. vtxFitter.setBeamSize(beamSize); vtxFitter.setBeamPosition(beamPositionToUse); @@ -880,11 +882,12 @@ private List shiftTracksToVertex(List parti double[] newRef = {vtxPos.z(), vtxPos.x(), 0.0};//the TrackUtils.getParametersAtNewRefPoint method only shifts in xy tracking frame List newTrks = new ArrayList(); for (ReconstructedParticle part : particles) { - BaseTrackState oldTS = (BaseTrackState) part.getTracks().get(0).getTrackStates().get(0); + BaseTrackState oldTS = (BaseTrackState) part.getTracks().get(0).getTrackStates().get(TrackState.AtPerigee); double[] newParams = TrackUtils.getParametersAtNewRefPoint(newRef, oldTS); SymmetricMatrix newCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, oldTS.getReferencePoint(), oldTS.getParameters(), new SymmetricMatrix(5, oldTS.getCovMatrix(), true)); //mg...I don't like this re-casting, but toBilliorTrack only takes Track as input - BaseTrackState newTS = new BaseTrackState(newParams, newRef, newCov.asPackedArray(true), TrackState.AtIP, bField); + //make the state with bfield = field at original z-position because TrackUtils.getParametersAtNewRefPoint does not change curvature! + BaseTrackState newTS = new BaseTrackState(newParams, newRef, newCov.asPackedArray(true), TrackState.AtPerigee, oldTS.getBLocal()); BilliorTrack electronBTrackShift = this.toBilliorTrack(newTS); newTrks.add(electronBTrackShift); }