diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000000..72074576ac --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,25 @@ +name: build + +on: + workflow_dispatch: + pull_request: + branches: [master] + +jobs: + build: + runs-on: ubuntu-latest + steps: + - name: Checkout source code + uses: actions/checkout@v3 + + - name: Set up JDK 21 + uses: actions/setup-java@v1 + with: + java-version: 21 + server-id: github + + - name: Build + run: | + mvn -B install -DskipTests=true + env: + GITHUB_TOKEN: ${{secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index cccecf95b7..90f87da06f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -12,11 +12,14 @@ jobs: - name: Checkout source code uses: actions/checkout@v3 - - name: Set up JDK 1.8 + - name: Set up JDK 21 uses: actions/setup-java@v1 with: - java-version: 1.8 + java-version: 21 + server-id: github - name: Build and run tests run: | mvn -B install -T4 + env: + GITHUB_TOKEN: ${{secrets.GITHUB_TOKEN }} diff --git a/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java b/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java index 95c223877f..bcd782d849 100644 --- a/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java +++ b/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java @@ -137,7 +137,7 @@ protected void process(EventHeader event) { truth_trk.getTrackStates().clear(); double[] ref = new double[] { 0., 0., 0. }; SymmetricMatrix cov = new SymmetricMatrix(5); - TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield); + TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield); truth_trk.getTrackStates().add(stateIP); truth_trk.setChisq(-1); truth_trk.setNDF(-1); 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/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java b/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java index b3615a4713..7673ff8ffa 100644 --- a/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java +++ b/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java @@ -803,8 +803,8 @@ else if (daug.getPDGID() == -11) double[] posShiftPars = TrackUtils.getParametersAtNewRefPoint(newRef, posOldTS); SymmetricMatrix eleShiftCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, eleOldTS.getReferencePoint(), eleOldTS.getParameters(), new SymmetricMatrix(5, eleOldTS.getCovMatrix(), true)); SymmetricMatrix posShiftCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, posOldTS.getReferencePoint(), posOldTS.getParameters(), new SymmetricMatrix(5, posOldTS.getCovMatrix(), true)); - BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleShiftCov.asPackedArray(true), TrackState.AtIP, B_FIELD); - BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posShiftCov.asPackedArray(true), TrackState.AtIP, B_FIELD); + BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleShiftCov.asPackedArray(true), TrackState.AtPerigee, B_FIELD); + BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posShiftCov.asPackedArray(true), TrackState.AtPerigee, B_FIELD); // BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleOldTS.getCovMatrix(), TrackState.AtIP, B_FIELD); // BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posOldTS.getCovMatrix(), TrackState.AtIP, B_FIELD); BilliorTrack electronBTrackShift = toBilliorTrack(eleShiftTS); diff --git a/integration-tests/src/test/java/org/hps/test/it/PhysRun2016FeeReconKalmanTest.java b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016FeeReconKalmanTest.java new file mode 100644 index 0000000000..d118c4f70a --- /dev/null +++ b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016FeeReconKalmanTest.java @@ -0,0 +1,216 @@ +package org.hps.test.it; + +import java.util.List; + +import org.hps.recon.ecal.cluster.ClusterUtilities; +import org.hps.recon.tracking.TrackData; +import org.hps.recon.tracking.TrackType; +import org.hps.recon.tracking.TrackUtils; +import org.hps.record.triggerbank.AbstractIntData; +import org.hps.record.triggerbank.TIData; +import org.lcsim.event.Cluster; +import org.lcsim.event.EventHeader; +import org.lcsim.event.GenericObject; +import org.lcsim.event.ReconstructedParticle; +import org.lcsim.event.Track; +import org.lcsim.math.chisq.ChisqProb; + +import hep.aida.IHistogram1D; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +/** + * Test FEE reconstruction for 2016 physics run + */ +public class PhysRun2016FeeReconKalmanTest extends ReconTest { + + static final String DETECTOR = "HPS-PhysicsRun2016-v5-3-fieldmap_v4_globalAlign"; + static final String TEST_FILE_NAME = "hps_007796_feeskim.evio"; + static final String STEERING = "/org/hps/steering/recon/PhysicsRun2016FullRecon_KF_TrackClusterMatcher.lcsim"; + static final int NEVENTS = 5000; + static final long MAX_EVENT_TIME = -1L; + + public PhysRun2016FeeReconKalmanTest() { + super(PhysRun2016FeeReconKalmanTest.class, + DETECTOR, + TEST_FILE_NAME, + STEERING, + NEVENTS, + new PlotDriver(), + DEFAULT_TOLERANCE, + MAX_EVENT_TIME, + true, + true); + } + + private static class PlotDriver extends RefDriver { + + IHistogram1D trkChisqNdfTop = aida.histogram1D("Top Track Chisq per DoF", 100, 0., 100.); + IHistogram1D trkChisqProbTop = aida.histogram1D("Top Track Chisq Prob", 100, 0., 1.); + IHistogram1D trkNhitsTop = aida.histogram1D("Top Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkMomentumTop = aida.histogram1D("Top Track Momentum", 200, 0., 3.); + IHistogram1D trkMomentumTopLT10 = aida.histogram1D("Top <10 Hit Track Momentum", 200, 0., 3.); + IHistogram1D trkMomentumTopGE10 = aida.histogram1D("Top >=10 Hit Track Momentum", 200, 0., 3.); + IHistogram1D trkdEdXTopLT10 = aida.histogram1D("Top <10 Hit Track dEdx", 100, 0., .00015); + IHistogram1D trkdEdXTopGE10 = aida.histogram1D("Top >=10 Hit Track dEdx", 100, 0., .00015); + IHistogram1D trkthetaTop = aida.histogram1D("Top Track theta", 100, 0.01, 0.05); + IHistogram1D trkX0Top = aida.histogram1D("Top Track X0", 100, -0.5, 0.5); + IHistogram1D trkY0Top = aida.histogram1D("Top Track Y0", 100, -5.0, 5.0); + IHistogram1D trkZ0Top = aida.histogram1D("Top Track Z0", 100, -1.0, 1.0); + + IHistogram1D trkChisqNdfBottom = aida.histogram1D("Bottom Track Chisq per DoF", 100, 0., 100.); + IHistogram1D trkChisqProbBottom = aida.histogram1D("Bottom Track Chisq Prob", 100, 0., 1.); + IHistogram1D trkNhitsBottom = aida.histogram1D("Bottom Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkMomentumBottom = aida.histogram1D("Bottom Track Momentum", 200, 0., 3.); + IHistogram1D trkMomentumBottomLT10 = aida.histogram1D("Bottom <10 Hit Track Momentum", 200, 0., 3.); + IHistogram1D trkMomentumBottomGE10 = aida.histogram1D("Bottom >=10 Hit Track Momentum", 200, 0., 3.); + IHistogram1D trkdEdXBottomLT10 = aida.histogram1D("Bottom <10 Hit Track dEdx", 100, 0., .00015); + IHistogram1D trkdEdXBottomGE10 = aida.histogram1D("Bottom >=10 Hit Track dEdx", 100, 0., .00015); + IHistogram1D trkthetaBottom = aida.histogram1D("Bottom Track theta", 100, 0.01, 0.05); + IHistogram1D trkX0Bottom = aida.histogram1D("Bottom Track X0", 100, -0.5, 0.5); + IHistogram1D trkY0Bottom = aida.histogram1D("Bottom Track Y0", 100, -5.0, 5.0); + IHistogram1D trkZ0Bottom = aida.histogram1D("Bottom Track Z0", 100, -1.0, 1.0); + + final String finalStateParticlesColName = "OtherElectrons"; + + Double beamEnergy = 2.306; + + // Set min seed energy value + final double seedCut = 1.2; + + // Set min cluster energy value + double clusterCut = 1.6; + + // Minimum number of hits per cluster + int minHits = 5; + + // Seed time cuts + double ctMin = 55.; + double ctMax = 61.; + + protected void process(EventHeader event) { + + super.process(event); + + // only keep singles triggers: + if (!event.hasCollection(GenericObject.class, "TriggerBank")) { + return; + } + boolean isSingles = false; + for (GenericObject gob : event.get(GenericObject.class, "TriggerBank")) { + if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) { + continue; + } + TIData tid = new TIData(gob); + if (tid.isSingle0Trigger() || tid.isSingle1Trigger()) { + isSingles = true; + break; + } + } + if (!isSingles) { + return; + } + if (!event.hasCollection(ReconstructedParticle.class, finalStateParticlesColName)) { + return; + } + List rpList = event.get(ReconstructedParticle.class, finalStateParticlesColName); + for (ReconstructedParticle rp : rpList) { + if (TrackType.isGBL(rp.getType())) { + continue; + } + if (rp.getMomentum().magnitude() > 1.5 * beamEnergy) { + continue; + } + // require both track and cluster + if (rp.getClusters().size() != 1) { + continue; + } + + if (rp.getTracks().size() != 1) { + continue; + } + + Track t = rp.getTracks().get(0); + double p = rp.getMomentum().magnitude(); + Cluster c = rp.getClusters().get(0); + // debug diagnostics to set cuts + if (debug) { + aida.cloud1D("clusterSeedHit energy").fill(ClusterUtilities.findSeedHit(c).getCorrectedEnergy()); + aida.cloud1D("cluster nHits").fill(c.getCalorimeterHits().size()); + aida.cloud2D("clusterSeedHit energy vs p").fill(p, ClusterUtilities.findSeedHit(c).getCorrectedEnergy()); + aida.cloud2D("cluster nHits vs p").fill(p, c.getCalorimeterHits().size()); + aida.cloud2D("cluster time vs p").fill(p, ClusterUtilities.getSeedHitTime(c)); + } + double ct = ClusterUtilities.getSeedHitTime(c); + + if (c.getEnergy() > clusterCut + && ClusterUtilities.findSeedHit(c).getCorrectedEnergy() > seedCut + && c.getCalorimeterHits().size() >= minHits + && ct > ctMin + && ct < ctMax) { + double chiSquared = t.getChi2(); + int ndf = t.getNDF(); + double chi2Ndf = t.getChi2() / t.getNDF(); + double chisqProb = ChisqProb.gammp(ndf, chiSquared); + int nHits = t.getTrackerHits().size(); + double dEdx = t.getdEdx(); + //rotate into physiscs frame of reference + Hep3Vector rprot = VecOp.mult(BEAM_AXIS_ROTATION, rp.getMomentum()); + double theta = Math.acos(rprot.z() / rprot.magnitude()); + + // debug diagnostics to set cuts + if (debug) { + aida.cloud1D("Track chisq per df").fill(chiSquared / ndf); + aida.cloud1D("Track chisq prob").fill(chisqProb); + //aida.cloud1D("Track nHits").fill(t.getTrackerHits().size()); + aida.cloud1D("Track momentum").fill(p); + aida.cloud1D("Track deDx").fill(t.getdEdx()); + aida.cloud1D("Track theta").fill(theta); + aida.cloud2D("Track theta vs p").fill(theta, p); + aida.cloud1D("rp x0").fill(TrackUtils.getX0(t)); + aida.cloud1D("rp y0").fill(TrackUtils.getY0(t)); + aida.cloud1D("rp z0").fill(TrackUtils.getZ0(t)); + } + + double trackDataTime = TrackData.getTrackTime(TrackData.getTrackData(event, t)); + aida.cloud1D("track data time").fill(trackDataTime); + if (isTopTrack(t)) { + trkChisqNdfTop.fill(chi2Ndf); + trkChisqProbTop.fill(chisqProb); + trkNhitsTop.fill(nHits); + trkMomentumTop.fill(p); + trkthetaTop.fill(theta); + trkX0Top.fill(TrackUtils.getX0(t)); + trkY0Top.fill(TrackUtils.getY0(t)); + trkZ0Top.fill(TrackUtils.getZ0(t)); + + if (nHits < 10) { + trkMomentumTopLT10.fill(p); + trkdEdXTopLT10.fill(dEdx); + } else if (nHits >=10 ) { + trkMomentumTopGE10.fill(p); + trkdEdXTopGE10.fill(dEdx); + } + } else { + trkChisqNdfBottom.fill(chi2Ndf); + trkChisqProbBottom.fill(chisqProb); + trkNhitsBottom.fill(nHits); + trkMomentumBottom.fill(p); + trkthetaBottom.fill(theta); + trkX0Bottom.fill(TrackUtils.getX0(t)); + trkY0Bottom.fill(TrackUtils.getY0(t)); + trkZ0Bottom.fill(TrackUtils.getZ0(t)); + + if (nHits<10 ) { + trkMomentumBottomLT10.fill(p); + trkdEdXBottomLT10.fill(dEdx); + } else if (nHits>=10) { + trkMomentumBottomGE10.fill(p); + trkdEdXBottomGE10.fill(dEdx); + } + } + } // end of cluster cuts + } + } + } +} diff --git a/integration-tests/src/test/java/org/hps/test/it/PhysRun2016MollerReconKalmanTest.java b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016MollerReconKalmanTest.java new file mode 100644 index 0000000000..f1ff2ffdbb --- /dev/null +++ b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016MollerReconKalmanTest.java @@ -0,0 +1,312 @@ +package org.hps.test.it; + +import static java.lang.Math.abs; + +import java.util.List; + +import org.hps.recon.tracking.TrackType; +import org.lcsim.event.EventHeader; +import org.lcsim.event.ReconstructedParticle; +import org.lcsim.event.Track; +import org.lcsim.event.Vertex; +import org.lcsim.math.chisq.ChisqProb; + +import hep.aida.IHistogram1D; +import hep.aida.IHistogram2D; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +/** + * Test Moller reconstruction for 2016 physics run + */ +public class PhysRun2016MollerReconKalmanTest extends ReconTest { + + static final String DETECTOR = "HPS-PhysicsRun2016-v5-3-fieldmap_v4_globalAlign"; + static final String TEST_FILE_NAME = "hps_007796_mollerskim.evio"; + static final String STEERING = "/org/hps/steering/recon/PhysicsRun2016FullRecon_KF_TrackClusterMatcher.lcsim"; + static final int NEVENTS = 5000; + static final long MAX_EVENT_TIME = -1L; + + public PhysRun2016MollerReconKalmanTest() { + super(PhysRun2016MollerReconKalmanTest.class, + DETECTOR, + TEST_FILE_NAME, + STEERING, + NEVENTS, + new PlotDriver(), + DEFAULT_TOLERANCE, + MAX_EVENT_TIME, + true, + true); + } + + private static class PlotDriver extends RefDriver { + + String[] vertexCollectionNames = { + "UnconstrainedMollerVertices", + "BeamspotConstrainedMollerVertices", + "TargetConstrainedMollerVertices"}; + + Double beamEnergy = 2.306; + double psumDelta = 0.06; + double thetaSumCut = 0.0475; + double trackChi2NdfCut = 8.; + + double psumMin = (1 - psumDelta) * beamEnergy; + double psumMax = (1 + psumDelta) * beamEnergy; + + IHistogram1D invMassHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Momentum", 200, 0., 3.0); + IHistogram1D pxHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_UnconstrainedMollerVertices = aida.histogram1D("UnconstrainedMollerVertices/Moller Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_UnconstrainedMollerVertices = aida.histogram2D("UnconstrainedMollerVertices/Moller p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_UnconstrainedMollerVertices = aida.histogram2D("UnconstrainedMollerVertices/Moller p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_UnconstrainedMollerVertices = aida.histogram2D("UnconstrainedMollerVertices/Moller theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_UnconstrainedMollerVertices = aida.histogram2D("UnconstrainedMollerVertices/Moller p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_UnconstrainedMollerVertices = aida.histogram2D("UnconstrainedMollerVertices/Moller vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + IHistogram1D invMassHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Momentum", 200, 0., 3.0); + IHistogram1D pxHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_BeamspotConstrainedMollerVertices = aida.histogram1D("BeamspotConstrainedMollerVertices/Moller Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_BeamspotConstrainedMollerVertices = aida.histogram2D("BeamspotConstrainedMollerVertices/Moller p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_BeamspotConstrainedMollerVertices = aida.histogram2D("BeamspotConstrainedMollerVertices/Moller p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_BeamspotConstrainedMollerVertices = aida.histogram2D("BeamspotConstrainedMollerVertices/Moller theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_BeamspotConstrainedMollerVertices = aida.histogram2D("BeamspotConstrainedMollerVertices/Moller p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_BeamspotConstrainedMollerVertices = aida.histogram2D("BeamspotConstrainedMollerVertices/Moller vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + IHistogram1D invMassHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Momentum", 200, 0., 3.0); + IHistogram1D pxHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_TargetConstrainedMollerVertices = aida.histogram1D("TargetConstrainedMollerVertices/Moller Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_TargetConstrainedMollerVertices = aida.histogram2D("TargetConstrainedMollerVertices/Moller p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_TargetConstrainedMollerVertices = aida.histogram2D("TargetConstrainedMollerVertices/Moller p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_TargetConstrainedMollerVertices = aida.histogram2D("TargetConstrainedMollerVertices/Moller theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_TargetConstrainedMollerVertices = aida.histogram2D("TargetConstrainedMollerVertices/Moller p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_TargetConstrainedMollerVertices = aida.histogram2D("TargetConstrainedMollerVertices/Moller vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + protected void process(EventHeader event) { + super.process(event); + for (String vertexCollectionName : vertexCollectionNames) { + + List vertices = event.get(Vertex.class, vertexCollectionName); + for (Vertex v : vertices) { + ReconstructedParticle rp = v.getAssociatedParticle(); + int type = rp.getType(); + boolean isGbl = TrackType.isGBL(type); + // require Kalman tracks in vertex + if (!isGbl) { + List parts = rp.getParticles(); + ReconstructedParticle rp1 = parts.get(0); + ReconstructedParticle rp2 = parts.get(1); + // basic sanity check here, remove full energy electrons (fee) + if (rp1.getMomentum().magnitude() > 1.5 * beamEnergy || rp2.getMomentum().magnitude() > 1.5 * beamEnergy) { + continue; + } + // require both reconstructed particles to have a track + // and at least one to have a cluster + if ((rp1.getClusters() == null) && (rp2.getClusters() == null)) { + continue; + } + + if (rp1.getTracks().size() != 1) { + continue; + } + if (rp2.getTracks().size() != 1) { + continue; + } + Track t1 = rp1.getTracks().get(0); + Track t2 = rp2.getTracks().get(0); + // require momentum sum to equal beam energy +- + double psum = rp1.getMomentum().magnitude() + rp2.getMomentum().magnitude(); + if (psum < psumMin || psum > psumMax) { + continue; + } + //rotate into physics frame of reference + Hep3Vector rprot = VecOp.mult(BEAM_AXIS_ROTATION, rp.getMomentum()); + Hep3Vector p1rot = VecOp.mult(BEAM_AXIS_ROTATION, rp1.getMomentum()); + Hep3Vector p2rot = VecOp.mult(BEAM_AXIS_ROTATION, rp2.getMomentum()); + double theta1 = Math.acos(p1rot.z() / p1rot.magnitude()); + double theta2 = Math.acos(p2rot.z() / p2rot.magnitude()); + double thetasum = theta1 + theta2; + // cut on thetasum + if (thetasum > thetaSumCut) { + continue; + } + // cut on Moller pX + if (abs(rprot.x()) > 0.01) { + continue; + } + // cut on Moller pY + if (abs(rp.getMomentum().y()) > .01) { + continue; + } + double t1ChisqNdf = t1.getChi2() / t1.getNDF(); + double t2ChisqNdf = t2.getChi2() / t2.getNDF(); + + double t1ChisqProb = ChisqProb.gammp(t1.getNDF(), t1.getChi2()); + double t2ChisqProb = ChisqProb.gammp(t2.getNDF(), t2.getChi2()); + // used to cut on prob < 0.995, which corresponds to roughly 3.4 + // change this to a cut on chi-squared/dof which people are more familiar with. + // Omar currently cuts on chi-squared <40(!), irrespective of 5 or 6 hit tracks + // let's start at chisq/dof of 8 + if (t1ChisqNdf > trackChi2NdfCut) {//(t1ChisqProb > 0.995) { + continue; + } + if (t2ChisqNdf > trackChi2NdfCut) {//(t2ChisqProb > 0.995) { + continue; + } + // all cuts passed, let's fill some histograms + Hep3Vector pos = v.getPosition(); + double p1 = rp1.getMomentum().magnitude(); + double p2 = rp2.getMomentum().magnitude(); + if (vertexCollectionName.equals("UnconstrainedMollerVertices")) { + invMassHist_UnconstrainedMollerVertices.fill(rp.getMass()); + pHist_UnconstrainedMollerVertices.fill(rp.getMomentum().magnitude()); + pxHist_UnconstrainedMollerVertices.fill(rprot.x()); + pyHist_UnconstrainedMollerVertices.fill(rprot.y()); + pzHist_UnconstrainedMollerVertices.fill(rprot.z()); + trkpHist_UnconstrainedMollerVertices.fill(p1); + trkpHist_UnconstrainedMollerVertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_UnconstrainedMollerVertices.fill(p1); + ptopvspbotHist_UnconstrainedMollerVertices.fill(p1, p2); + } else { + trkpbotHist_UnconstrainedMollerVertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_UnconstrainedMollerVertices.fill(p2); + } else { + trkpbotHist_UnconstrainedMollerVertices.fill(p2); + } + trkNhitsHist_UnconstrainedMollerVertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_UnconstrainedMollerVertices.fill(t2.getTrackerHits().size()); + trkChisqHist_UnconstrainedMollerVertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_UnconstrainedMollerVertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_UnconstrainedMollerVertices.fill(t1ChisqProb); + trkChisqProbHist_UnconstrainedMollerVertices.fill(t2ChisqProb); + vtxXHist_UnconstrainedMollerVertices.fill(pos.x()); + vtxYHist_UnconstrainedMollerVertices.fill(pos.y()); + vtxZHist_UnconstrainedMollerVertices.fill(pos.z()); + vtxChisqHist_UnconstrainedMollerVertices.fill(v.getChi2()); + + p1vsp2Hist_UnconstrainedMollerVertices.fill(p1, p2); + theta1vstheta2Hist_UnconstrainedMollerVertices.fill(theta1, theta2); + pvsthetaHist_UnconstrainedMollerVertices.fill(p1, theta1); + pvsthetaHist_UnconstrainedMollerVertices.fill(p2, theta2); + xvsyHist_UnconstrainedMollerVertices.fill(pos.x(), pos.y()); + } + if (vertexCollectionName.equals("BeamspotConstrainedMollerVertices")) { + invMassHist_BeamspotConstrainedMollerVertices.fill(rp.getMass()); + pHist_BeamspotConstrainedMollerVertices.fill(rp.getMomentum().magnitude()); + pxHist_BeamspotConstrainedMollerVertices.fill(rprot.x()); + pyHist_BeamspotConstrainedMollerVertices.fill(rprot.y()); + pzHist_BeamspotConstrainedMollerVertices.fill(rprot.z()); + trkpHist_BeamspotConstrainedMollerVertices.fill(p1); + trkpHist_BeamspotConstrainedMollerVertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_BeamspotConstrainedMollerVertices.fill(p1); + ptopvspbotHist_BeamspotConstrainedMollerVertices.fill(p1, p2); + } else { + trkpbotHist_BeamspotConstrainedMollerVertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_BeamspotConstrainedMollerVertices.fill(p2); + } else { + trkpbotHist_BeamspotConstrainedMollerVertices.fill(p2); + } + trkNhitsHist_BeamspotConstrainedMollerVertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_BeamspotConstrainedMollerVertices.fill(t2.getTrackerHits().size()); + trkChisqHist_BeamspotConstrainedMollerVertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_BeamspotConstrainedMollerVertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_BeamspotConstrainedMollerVertices.fill(t1ChisqProb); + trkChisqProbHist_BeamspotConstrainedMollerVertices.fill(t2ChisqProb); + vtxXHist_BeamspotConstrainedMollerVertices.fill(pos.x()); + vtxYHist_BeamspotConstrainedMollerVertices.fill(pos.y()); + vtxZHist_BeamspotConstrainedMollerVertices.fill(pos.z()); + vtxChisqHist_BeamspotConstrainedMollerVertices.fill(v.getChi2()); + + p1vsp2Hist_BeamspotConstrainedMollerVertices.fill(p1, p2); + theta1vstheta2Hist_BeamspotConstrainedMollerVertices.fill(theta1, theta2); + pvsthetaHist_BeamspotConstrainedMollerVertices.fill(p1, theta1); + pvsthetaHist_BeamspotConstrainedMollerVertices.fill(p2, theta2); + xvsyHist_BeamspotConstrainedMollerVertices.fill(pos.x(), pos.y()); + } + if (vertexCollectionName.equals("TargetConstrainedMollerVertices")) { + invMassHist_TargetConstrainedMollerVertices.fill(rp.getMass()); + pHist_TargetConstrainedMollerVertices.fill(rp.getMomentum().magnitude()); + pxHist_TargetConstrainedMollerVertices.fill(rprot.x()); + pyHist_TargetConstrainedMollerVertices.fill(rprot.y()); + pzHist_TargetConstrainedMollerVertices.fill(rprot.z()); + trkpHist_TargetConstrainedMollerVertices.fill(p1); + trkpHist_TargetConstrainedMollerVertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_TargetConstrainedMollerVertices.fill(p1); + ptopvspbotHist_TargetConstrainedMollerVertices.fill(p1, p2); + } else { + trkpbotHist_TargetConstrainedMollerVertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_TargetConstrainedMollerVertices.fill(p2); + } else { + trkpbotHist_TargetConstrainedMollerVertices.fill(p2); + } + trkNhitsHist_TargetConstrainedMollerVertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_TargetConstrainedMollerVertices.fill(t2.getTrackerHits().size()); + trkChisqHist_TargetConstrainedMollerVertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_TargetConstrainedMollerVertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_TargetConstrainedMollerVertices.fill(t1ChisqProb); + trkChisqProbHist_TargetConstrainedMollerVertices.fill(t2ChisqProb); + vtxXHist_TargetConstrainedMollerVertices.fill(pos.x()); + vtxYHist_TargetConstrainedMollerVertices.fill(pos.y()); + vtxZHist_TargetConstrainedMollerVertices.fill(pos.z()); + vtxChisqHist_TargetConstrainedMollerVertices.fill(v.getChi2()); + + p1vsp2Hist_TargetConstrainedMollerVertices.fill(p1, p2); + theta1vstheta2Hist_TargetConstrainedMollerVertices.fill(theta1, theta2); + pvsthetaHist_TargetConstrainedMollerVertices.fill(p1, theta1); + pvsthetaHist_TargetConstrainedMollerVertices.fill(p2, theta2); + xvsyHist_TargetConstrainedMollerVertices.fill(pos.x(), pos.y()); + } + } //Loop over Kalman-vertices + } // Loop over vertices + } // loop over various vertex collections + } + } +} diff --git a/integration-tests/src/test/java/org/hps/test/it/PhysRun2016V0ReconKalmanTest.java b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016V0ReconKalmanTest.java new file mode 100644 index 0000000000..c739cf4206 --- /dev/null +++ b/integration-tests/src/test/java/org/hps/test/it/PhysRun2016V0ReconKalmanTest.java @@ -0,0 +1,300 @@ +package org.hps.test.it; + +import static java.lang.Math.abs; + +import java.util.List; + +import org.hps.recon.ecal.cluster.ClusterUtilities; +import org.hps.recon.tracking.TrackType; +import org.lcsim.event.Cluster; +import org.lcsim.event.EventHeader; +import org.lcsim.event.ReconstructedParticle; +import org.lcsim.event.Track; +import org.lcsim.event.Vertex; +import org.lcsim.math.chisq.ChisqProb; + +import hep.aida.IHistogram1D; +import hep.aida.IHistogram2D; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +/** + * Test V0 reconstruction for 2016 physics run + */ +public class PhysRun2016V0ReconKalmanTest extends ReconTest { + + static final String DETECTOR = "HPS-PhysicsRun2016-v5-3-fieldmap_v4_globalAlign"; + static final String TEST_FILE_NAME = "hps_007796_v0skim.evio"; + static final String STEERING = "/org/hps/steering/recon/PhysicsRun2016FullRecon_KF_TrackClusterMatcher.lcsim"; + static final int NEVENTS = 5000; + static final long MAX_EVENT_TIME = -1L; + + public PhysRun2016V0ReconKalmanTest() { + super(PhysRun2016V0ReconKalmanTest.class, + DETECTOR, + TEST_FILE_NAME, + STEERING, + NEVENTS, + new PlotDriver(), + DEFAULT_TOLERANCE, + MAX_EVENT_TIME, + true, + true); + } + + private static class PlotDriver extends RefDriver { + + String[] vertexCollectionNames = { + "UnconstrainedV0Vertices_KF", "BeamspotConstrainedV0Vertices_KF", "TargetConstrainedV0Vertices_KF"}; + + Double beamEnergy = 2.306; + double trackChi2NdfCut = 100.; //corresponds to chisquared cut of 40 for 5-hit tracks + + IHistogram1D invMassHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Momentum", 200, 0., 3.0); + IHistogram1D pxHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxZHistL1L1_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Vertex z L1L1", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_UnconstrainedV0Vertices = aida.histogram1D("UnconstrainedV0Vertices/V0 Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_UnconstrainedV0Vertices = aida.histogram2D("UnconstrainedV0Vertices/V0 p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_UnconstrainedV0Vertices = aida.histogram2D("UnconstrainedV0Vertices/V0 p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_UnconstrainedV0Vertices = aida.histogram2D("UnconstrainedV0Vertices/V0 theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_UnconstrainedV0Vertices = aida.histogram2D("UnconstrainedV0Vertices/V0 p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_UnconstrainedV0Vertices = aida.histogram2D("UnconstrainedV0Vertices/V0 vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + IHistogram1D invMassHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Momentum", 200, 0., 3.0); + IHistogram1D pxHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_BeamspotConstrainedV0Vertices = aida.histogram1D("BeamspotConstrainedV0Vertices/V0 Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_BeamspotConstrainedV0Vertices = aida.histogram2D("BeamspotConstrainedV0Vertices/V0 p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_BeamspotConstrainedV0Vertices = aida.histogram2D("BeamspotConstrainedV0Vertices/V0 p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_BeamspotConstrainedV0Vertices = aida.histogram2D("BeamspotConstrainedV0Vertices/V0 theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_BeamspotConstrainedV0Vertices = aida.histogram2D("BeamspotConstrainedV0Vertices/V0 p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_BeamspotConstrainedV0Vertices = aida.histogram2D("BeamspotConstrainedV0Vertices/V0 vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + IHistogram1D invMassHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Invariant Mass", 200, 0., 0.1); + IHistogram1D pHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Momentum", 200, 0., 3.0); + IHistogram1D pxHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 x Momentum", 200, -0.01, 0.01); + IHistogram1D pyHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 y Momentum", 200, -0.01, 0.01); + IHistogram1D pzHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 z Momentum", 200, 0., 3.0); + IHistogram1D trkpHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkptopHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Top Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkpbotHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Bottom Track Momentum", 100, 0.25, 1.75); + IHistogram1D trkNhitsHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Track Number of Hits", 15, -0.5, 14.5); + IHistogram1D trkChisqHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Track Chisq per DoF", 100, 0.0, 20.0); + IHistogram1D trkChisqProbHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Track Chisq Prob", 100, 0.0, 1.0); + IHistogram1D vtxXHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Vertex x", 200, -2.5, 2.5); + IHistogram1D vtxYHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Vertex y", 200, -1.0, 1.0); + IHistogram1D vtxZHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Vertex z", 200, -20.0, 20.0); + IHistogram1D vtxChisqHist_TargetConstrainedV0Vertices = aida.histogram1D("TargetConstrainedV0Vertices/V0 Vertex Chisq", 100, 0.0, 100.0); + + IHistogram2D p1vsp2Hist_TargetConstrainedV0Vertices = aida.histogram2D("TargetConstrainedV0Vertices/V0 p1 vs p2", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D ptopvspbotHist_TargetConstrainedV0Vertices = aida.histogram2D("TargetConstrainedV0Vertices/V0 p top vs p bottom", 200, 0.25, 1.75, 200, 0.25, 1.75); + IHistogram2D theta1vstheta2Hist_TargetConstrainedV0Vertices = aida.histogram2D("TargetConstrainedV0Vertices/V0 theta1 vs theta2", 100, 0.01, 0.05, 100, 0.01, 0.05); + IHistogram2D pvsthetaHist_TargetConstrainedV0Vertices = aida.histogram2D("TargetConstrainedV0Vertices/V0 p vs theta", 100, 0.25, 1.75, 100, 0.01, 0.05); + IHistogram2D xvsyHist_TargetConstrainedV0Vertices = aida.histogram2D("TargetConstrainedV0Vertices/V0 vertex X vs Y", 250, -2.5, 2.5, 100, -1.0, 1.0); + + protected void process(EventHeader event) { + super.process(event); + for (String vertexCollectionName : vertexCollectionNames) { + List vertices = event.get(Vertex.class, vertexCollectionName); + for (Vertex v : vertices) { + ReconstructedParticle rp = v.getAssociatedParticle(); + int type = rp.getType(); + boolean isGbl = TrackType.isGBL(type); + // require Kalman tracks in vertex + if (!isGbl) { + List parts = rp.getParticles(); + ReconstructedParticle rp1 = parts.get(0); + ReconstructedParticle rp2 = parts.get(1); + // basic sanity check here, remove full energy electrons (fee) + if (rp1.getMomentum().magnitude() > 1.5 * beamEnergy || rp2.getMomentum().magnitude() > 1.5 * beamEnergy) { + continue; + } + // require both reconstructed particles to have a track and a cluster + if (rp1.getClusters().size() != 1) { + continue; + } + if (rp2.getClusters().size() != 1) { + continue; + } + if (rp1.getTracks().size() != 1) { + continue; + } + if (rp2.getTracks().size() != 1) { + continue; + } + Track t1 = rp1.getTracks().get(0); + Track t2 = rp2.getTracks().get(0); + Cluster c1 = rp1.getClusters().get(0); + Cluster c2 = rp2.getClusters().get(0); + double deltaT = ClusterUtilities.getSeedHitTime(c1) - ClusterUtilities.getSeedHitTime(c2); + // require cluster times to be coincident within 2 ns + if (abs(deltaT) > 2.0) { + continue; + } + //rotate into physiscs frame of reference + Hep3Vector rprot = VecOp.mult(BEAM_AXIS_ROTATION, rp.getMomentum()); + Hep3Vector p1rot = VecOp.mult(BEAM_AXIS_ROTATION, rp1.getMomentum()); + Hep3Vector p2rot = VecOp.mult(BEAM_AXIS_ROTATION, rp2.getMomentum()); + double theta1 = Math.acos(p1rot.z() / p1rot.magnitude()); + double theta2 = Math.acos(p2rot.z() / p2rot.magnitude()); + double t1ChisqNdf = t1.getChi2() / t1.getNDF(); + double t2ChisqNdf = t2.getChi2() / t2.getNDF(); + + double t1ChisqProb = ChisqProb.gammp(t1.getNDF(), t1.getChi2()); + double t2ChisqProb = ChisqProb.gammp(t2.getNDF(), t2.getChi2()); + // used to cut on prob < 0.995, which corresponds to roughly 3.4 + // change this to a cut on chi-squared/dof which people are more familiar with. + // Omar currently cuts on chi-squared <40(!), irrespective of 5 or 6 hit tracks + // let's start at chisq/dof of 8 + if (t1ChisqNdf > trackChi2NdfCut) {//(t1ChisqProb > 0.995) { + continue; + } + if (t2ChisqNdf > trackChi2NdfCut) {//(t2ChisqProb > 0.995) { + continue; + } + // all cuts passed, let's fill some histograms + Hep3Vector pos = v.getPosition(); + double p1 = rp1.getMomentum().magnitude(); + double p2 = rp2.getMomentum().magnitude(); + if (vertexCollectionName.equals("UnconstrainedV0Vertices_KF")) { + invMassHist_UnconstrainedV0Vertices.fill(rp.getMass()); + pHist_UnconstrainedV0Vertices.fill(rp.getMomentum().magnitude()); + pxHist_UnconstrainedV0Vertices.fill(rprot.x()); + pyHist_UnconstrainedV0Vertices.fill(rprot.y()); + pzHist_UnconstrainedV0Vertices.fill(rprot.z()); + trkpHist_UnconstrainedV0Vertices.fill(p1); + trkpHist_UnconstrainedV0Vertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_UnconstrainedV0Vertices.fill(p1); + ptopvspbotHist_UnconstrainedV0Vertices.fill(p1, p2); + } else { + trkpbotHist_UnconstrainedV0Vertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_UnconstrainedV0Vertices.fill(p2); + } else { + trkpbotHist_UnconstrainedV0Vertices.fill(p2); + } + trkNhitsHist_UnconstrainedV0Vertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_UnconstrainedV0Vertices.fill(t2.getTrackerHits().size()); + trkChisqHist_UnconstrainedV0Vertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_UnconstrainedV0Vertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_UnconstrainedV0Vertices.fill(t1ChisqProb); + trkChisqProbHist_UnconstrainedV0Vertices.fill(t2ChisqProb); + vtxXHist_UnconstrainedV0Vertices.fill(pos.x()); + vtxYHist_UnconstrainedV0Vertices.fill(pos.y()); + vtxZHist_UnconstrainedV0Vertices.fill(pos.z()); + if (hasLayer1Hit(t1) && hasLayer1Hit(t2)) { + vtxZHistL1L1_UnconstrainedV0Vertices.fill(pos.z()); + } + vtxChisqHist_UnconstrainedV0Vertices.fill(v.getChi2()); + + p1vsp2Hist_UnconstrainedV0Vertices.fill(p1, p2); + theta1vstheta2Hist_UnconstrainedV0Vertices.fill(theta1, theta2); + pvsthetaHist_UnconstrainedV0Vertices.fill(p1, theta1); + pvsthetaHist_UnconstrainedV0Vertices.fill(p2, theta2); + xvsyHist_UnconstrainedV0Vertices.fill(pos.x(), pos.y()); + } + if (vertexCollectionName.equals("BeamspotConstrainedV0Vertices_KF")) { + invMassHist_BeamspotConstrainedV0Vertices.fill(rp.getMass()); + pHist_BeamspotConstrainedV0Vertices.fill(rp.getMomentum().magnitude()); + pxHist_BeamspotConstrainedV0Vertices.fill(rprot.x()); + pyHist_BeamspotConstrainedV0Vertices.fill(rprot.y()); + pzHist_BeamspotConstrainedV0Vertices.fill(rprot.z()); + trkpHist_BeamspotConstrainedV0Vertices.fill(p1); + trkpHist_BeamspotConstrainedV0Vertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_BeamspotConstrainedV0Vertices.fill(p1); + ptopvspbotHist_BeamspotConstrainedV0Vertices.fill(p1, p2); + } else { + trkpbotHist_BeamspotConstrainedV0Vertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_BeamspotConstrainedV0Vertices.fill(p2); + } else { + trkpbotHist_BeamspotConstrainedV0Vertices.fill(p2); + } + trkNhitsHist_BeamspotConstrainedV0Vertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_BeamspotConstrainedV0Vertices.fill(t2.getTrackerHits().size()); + trkChisqHist_BeamspotConstrainedV0Vertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_BeamspotConstrainedV0Vertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_BeamspotConstrainedV0Vertices.fill(t1ChisqProb); + trkChisqProbHist_BeamspotConstrainedV0Vertices.fill(t2ChisqProb); + vtxXHist_BeamspotConstrainedV0Vertices.fill(pos.x()); + vtxYHist_BeamspotConstrainedV0Vertices.fill(pos.y()); + vtxZHist_BeamspotConstrainedV0Vertices.fill(pos.z()); + vtxChisqHist_BeamspotConstrainedV0Vertices.fill(v.getChi2()); + + p1vsp2Hist_BeamspotConstrainedV0Vertices.fill(p1, p2); + theta1vstheta2Hist_BeamspotConstrainedV0Vertices.fill(theta1, theta2); + pvsthetaHist_BeamspotConstrainedV0Vertices.fill(p1, theta1); + pvsthetaHist_BeamspotConstrainedV0Vertices.fill(p2, theta2); + xvsyHist_BeamspotConstrainedV0Vertices.fill(pos.x(), pos.y()); + } + if (vertexCollectionName.equals("TargetConstrainedV0Vertices_KF")) { + invMassHist_TargetConstrainedV0Vertices.fill(rp.getMass()); + pHist_TargetConstrainedV0Vertices.fill(rp.getMomentum().magnitude()); + pxHist_TargetConstrainedV0Vertices.fill(rprot.x()); + pyHist_TargetConstrainedV0Vertices.fill(rprot.y()); + pzHist_TargetConstrainedV0Vertices.fill(rprot.z()); + trkpHist_TargetConstrainedV0Vertices.fill(p1); + trkpHist_TargetConstrainedV0Vertices.fill(p2); + if (isTopTrack(t1)) { + trkptopHist_TargetConstrainedV0Vertices.fill(p1); + ptopvspbotHist_TargetConstrainedV0Vertices.fill(p1, p2); + } else { + trkpbotHist_TargetConstrainedV0Vertices.fill(p1); + } + if (isTopTrack(t2)) { + trkptopHist_TargetConstrainedV0Vertices.fill(p2); + } else { + trkpbotHist_TargetConstrainedV0Vertices.fill(p2); + } + trkNhitsHist_TargetConstrainedV0Vertices.fill(t1.getTrackerHits().size()); + trkNhitsHist_TargetConstrainedV0Vertices.fill(t2.getTrackerHits().size()); + trkChisqHist_TargetConstrainedV0Vertices.fill(t1.getChi2() / t1.getNDF()); + trkChisqHist_TargetConstrainedV0Vertices.fill(t2.getChi2() / t2.getNDF()); + trkChisqProbHist_TargetConstrainedV0Vertices.fill(t1ChisqProb); + trkChisqProbHist_TargetConstrainedV0Vertices.fill(t2ChisqProb); + vtxXHist_TargetConstrainedV0Vertices.fill(pos.x()); + vtxYHist_TargetConstrainedV0Vertices.fill(pos.y()); + vtxZHist_TargetConstrainedV0Vertices.fill(pos.z()); + vtxChisqHist_TargetConstrainedV0Vertices.fill(v.getChi2()); + + p1vsp2Hist_TargetConstrainedV0Vertices.fill(p1, p2); + theta1vstheta2Hist_TargetConstrainedV0Vertices.fill(theta1, theta2); + pvsthetaHist_TargetConstrainedV0Vertices.fill(p1, theta1); + pvsthetaHist_TargetConstrainedV0Vertices.fill(p2, theta2); + xvsyHist_TargetConstrainedV0Vertices.fill(pos.x(), pos.y()); + } + } //Loop over Kalman-vertices + } // Loop over vertices + } // loop over various vertex collections + } + } +} diff --git a/pom.xml b/pom.xml index 8f4bc25053..b79670e95d 100644 --- a/pom.xml +++ b/pom.xml @@ -19,7 +19,7 @@ false true true - 4.4.0 + 4.5.0 14.1 4.4.6 4.13.1 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..b10c5170f7 100644 --- a/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java +++ b/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java @@ -16,6 +16,7 @@ import org.hps.recon.vertexing.BilliorVertex; import org.hps.recon.vertexing.BilliorVertexer; import org.hps.record.StandardCuts; +import org.hps.recon.tracking.TrackStateUtils; import org.lcsim.detector.tracker.silicon.HpsSiSensor; import org.lcsim.event.EventHeader; import org.lcsim.event.RawTrackerHit; @@ -509,12 +510,22 @@ 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(); + double bLocal=TrackStateUtils.getTrackStatesAtLocation(electron.getTracks().get(0),TrackState.AtPerigee).get(0).getBLocal(); + // System.out.println("For Vertexer: "+TrackStateUtils.getTrackStatesAtLocation(electron.getTracks().get(0),TrackState.AtPerigee).get(0).toString()); + // if we are using GBL tracks (trackType=0), set bLocal to bField (i.e. at SVT center) + if(trackType==0) + bLocal=bField; + + //set up vertexer with field + // System.out.println("track type = "+trackType+"; setting vertex field = "+bLocal); + BilliorVertexer vtxFitter = new BilliorVertexer(bLocal); // TODO: The beam size should come from the conditions database. vtxFitter.setBeamSize(beamSize); vtxFitter.setBeamPosition(beamPositionToUse); @@ -880,11 +891,13 @@ 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)TrackStateUtils.getTrackStatesAtLocation(part.getTracks().get(0),TrackState.AtPerigee).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); } diff --git a/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java b/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java index 8d7aaf3811..a6a798c89d 100644 --- a/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java +++ b/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java @@ -153,6 +153,10 @@ public void setApplyClusterCorrections(boolean val) { * LCIO collection name for tracks. */ protected String matcherTrackCollectionName = "GBLTracks"; + /** + * track type: Kalman = 1; GBL = 0 obtained from matcherTrackCollectionName + */ + protected int trackType = 0; /** * Track Cluster Algorithm set to Kalman or GBL Tracks */ @@ -446,6 +450,10 @@ protected void detectorChanged(Detector detector) { matcher.setTrackCollectionName(matcherTrackCollectionName); matcher.enablePlots(enableTrackClusterMatchPlots); + //set the track type; default is 0 (GBL) + if(matcherTrackCollectionName.contains("Kalman") || matcherTrackCollectionName.contains("KF")) + trackType=1; + // Set the magnetic field parameters to the appropriate values. Hep3Vector ip = new BasicHep3Vector(0., 0., 500.0); bField = detector.getFieldMap().getField(ip).y(); diff --git a/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java b/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java index 4d4ebe745c..b39e129f12 100644 --- a/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java +++ b/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java @@ -1479,7 +1479,7 @@ protected void process(EventHeader event) { truth_trk.getTrackStates().clear(); double[] ref = new double[] { 0., 0., 0. }; SymmetricMatrix cov = new SymmetricMatrix(5); - TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield); + TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield); truth_trk.getTrackStates().add(stateIP); truth_trk.setChisq(-1); truth_trk.setNDF(-1); @@ -1538,7 +1538,7 @@ public double[] getMCPTrackParameters(MCParticle mcp, double bfield){ truth_trk.getTrackStates().clear(); double[] ref = new double[] { 0., 0., 0. }; SymmetricMatrix cov = new SymmetricMatrix(5); - TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield); + TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield); truth_trk.getTrackStates().add(stateIP); truth_trk.setChisq(-1); truth_trk.setNDF(-1); diff --git a/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java b/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java index 048c1d1a2b..aa2e33f96a 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java @@ -214,9 +214,9 @@ public void process(EventHeader event) { // Extrapolate the tweaked track to the face of the Ecal and get the // track state - TrackState stateIP = TrackUtils.getTrackStateAtLocation(track, TrackState.AtIP); + TrackState stateIP = TrackUtils.getTrackStateAtLocation(track, TrackState.AtPerigee); if (stateIP == null) { - throw new RuntimeException("IP track state for GBL track was not found"); + throw new RuntimeException("Perigee track state for GBL track was not found"); } TrackState stateEcalIP = TrackUtils.extrapolateTrackUsingFieldMap(stateIP, extStartPos, ecalPosition, stepSize, bFieldMap); diff --git a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java index 71ef5d4b03..bd0b9f3cd2 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java +++ b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java @@ -1750,18 +1750,18 @@ public static Double[] getIsolations(Track trk, RelationalTable hitToStrips, Rel * . */ public static BaseTrackState extrapolateTrackUsingFieldMap(Track track, double startPositionX, double endPositionX, double stepSize, FieldMap fieldMap) { - TrackState stateAtIP = null; + TrackState stateAtPerigee = null; for (TrackState state : track.getTrackStates()) { - if (state.getLocation() == TrackState.AtIP) { - stateAtIP = state; + if (state.getLocation() == TrackState.AtPerigee) { + stateAtPerigee = state; } } - if (stateAtIP == null) { - throw new RuntimeException("No track state at IP was found so this function shouldn't be used."); + if (stateAtPerigee == null) { + throw new RuntimeException("No track state at Perigee was found so this function shouldn't be used."); } // Extrapolate this track state - return extrapolateTrackUsingFieldMap(stateAtIP, startPositionX, endPositionX, stepSize, fieldMap); + return extrapolateTrackUsingFieldMap(stateAtPerigee, startPositionX, endPositionX, stepSize, fieldMap); } /** diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java index c198dbed8f..82b0092666 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java +++ b/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java @@ -123,9 +123,14 @@ public static Pair makeCorrectedTrack(FittedGblTrajectory fi trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge trk.getTrackStates().clear(); - TrackState stateIP = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield); - trk.getTrackStates().add(stateIP); + // Set the track state to AtPerigee to agree with new track state location definitions as of hps-lcsim version 4.5.0 + TrackState statePerigee = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtPerigee, bfield); + trk.getTrackStates().add(statePerigee); + //also add state at the IP for backward compatibility, but it's the same as the perigee for gbl tracks + TrackState stateIP = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield); + trk.getTrackStates().add(stateIP); + if (!storeTrackStates) { // just store last state Pair correctedHelixParamsLast = fittedGblTrajectory.getCorrectedPerigeeParameters(helicalTrackFit, FittedGblTrajectory.GBLPOINT.LAST, bfield); 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..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 @@ -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 @@ -393,6 +395,7 @@ HelixState propagateRungeKutta(Plane pln, 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/KFOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java index c1a7b486a8..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,10 @@ 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; @@ -263,10 +264,15 @@ public void process(EventHeader event) { particles = event.get(ReconstructedParticle.class, inputCollectionName); for (ReconstructedParticle particle : particles) { //this requires track cluster match + if(debug) + System.out.println(this.getClass().getName()+":: from ReconParticle found "+particle.getTracks().size()+" tracks and "+ + particle.getClusters().size()+" clusters "); 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); } @@ -280,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; @@ -290,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; @@ -333,13 +343,10 @@ 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); - // if (b_doGBLkinks) - // doGBLkinks(trk,gblKink, sensorNums); - if (useParticles) doEoPPlots(trk,TrackClusterPairs.get(trk)); @@ -372,15 +379,15 @@ private void doEoPPlots(Track track, Cluster cluster) { double trkCluTime=trackTime-cluster.getCalorimeterHits().get(0).getTime()-timeOffset; aidaKF.histogram1D(trkpFolder+"trk-cluTime_"+charge+"_"+vol).fill(trkCluTime); aidaKF.histogram1D(trkpFolder+"trk-cluTime_"+vol).fill(trkCluTime); - - - aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol).fill(trackp,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol).fill(tanL,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol).fill(phi,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(tanL,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_phi").fill(phi,eop); - // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol).fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol).fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol).fill(phi,eop); + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi").fill(phi,eop); + + // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, // phi, // eop); @@ -389,17 +396,19 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaKF.histogram1D(eopFolder+"Ecluster_"+vol+"_fid").fill(energy); aidaKF.histogram1D(eopFolder+"EoP_"+vol+"_fid").fill(eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+vol+"_fid").fill(phi,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+vol+"_fid").fill(trackp,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+vol+"_fid").fill(tanL,eop); + if(b_doBig2DPlots){ + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+vol+"_fid").fill(phi,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+vol+"_fid").fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+vol+"_fid").fill(tanL,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol+"_fid").fill(trackp,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol+"_fid").fill(tanL,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol+"_fid").fill(phi,eop); - - aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(tanL,eop); - aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(phi,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol+"_fid").fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol+"_fid").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol+"_fid").fill(phi,eop); + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(phi,eop); + } // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(tanL, //phi, //eop); @@ -425,36 +434,32 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaKF.histogram1D(eopFolder+"trk_clu_resX_"+vol+"_fid").fill(trkX-clusterX); aidaKF.histogram1D(eopFolder+"trk_clu_resY_"+vol+"_fid").fill(trkY-clusterY); - aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX_"+vol+"_fid").fill(trkX,trkX-clusterX); aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY_"+vol+"_fid").fill(trkY,trkX-clusterX); aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX_"+vol+"_fid").fill(trkX,trkY-clusterY); aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY_"+vol+"_fid").fill(trkY,trkY-clusterY); - + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+vol+"_fid").fill(trackp,trkY-clusterY); aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+vol+"_fid").fill(trackp,trkX-clusterX); aidaKF.histogram2D(eopFolder+"trkY_vs_tanL_"+vol+"_fid").fill(tanL,trkY); - aidaKF.histogram1D(eopFolder+"Xcluster_"+charge+"_"+vol+"_fid").fill(clusterX); aidaKF.histogram1D(eopFolder+"Ycluster_"+charge+"_"+vol+"_fid").fill(clusterY); aidaKF.histogram1D(eopFolder+"trk_clu_resX_"+charge+"_"+vol+"_fid").fill(trkX-clusterX); aidaKF.histogram1D(eopFolder+"trk_clu_resY_"+charge+"_"+vol+"_fid").fill(trkY-clusterY); - aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkX-clusterX); aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkX-clusterX); aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkY-clusterY); aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkY-clusterY); - + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkY-clusterY); aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkX-clusterX); - + aidaKF.histogram2D(eopFolder+"trkY_vs_tanL_"+charge+"_"+vol+"_fid").fill(tanL,trkY); - // @@ -530,8 +535,10 @@ private void FillKFTrackPlot(String str, String isTop, String charge, double val } private void FillKFTrackPlot(String str, String isTop, String charge, double valX, double valY) { - aidaKF.histogram2D(str+isTop).fill(valX,valY); - aidaKF.histogram2D(str+isTop+charge).fill(valX,valY); + if(b_doBig2DPlots){ + aidaKF.histogram2D(str+isTop).fill(valX,valY); + aidaKF.histogram2D(str+isTop+charge).fill(valX,valY); + } } /* @@ -662,21 +669,22 @@ private void doBasicKFtrack(Track trk, Map 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 +754,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 +770,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 +903,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 +920,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 +989,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 +1088,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 +1131,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 +1233,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 74998a3a40..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 @@ -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 */ @@ -710,6 +710,39 @@ public double originArcLength() { } return arcLength[0]; } + + public HelixState getHelixAtPlane(Plane xPlane, boolean pivotAtIntersect){ + 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..e72b67dc69 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,8 +107,9 @@ 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; + private boolean saveTrackStateAtIntercept = true; public void setRunNumber(int runNumber){ this.runNumber = runNumber; @@ -254,7 +255,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(); } @@ -311,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; @@ -344,7 +355,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. @@ -356,14 +368,28 @@ 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); + 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) { + 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) { + // 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)) { @@ -405,38 +431,47 @@ public TrackState createTrackState(MeasurementSite ms, int loc, boolean useSmoot // 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) { + // 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 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, posTrking, covHPS, loc, bLocal); } public void printGBLStripClusterData(GBLStripClusterData clstr) { @@ -571,23 +606,24 @@ 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(); - // 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 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[] 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, alpha); + double[] newCov = getLCSimCov(globalCov, alpha).asPackedArray(true); + TrackState ts = new BaseTrackState(newParams, refD ,newCov, TrackState.AtPerigee, B); if (ts != null) { newTrack.getTrackStates().add(ts); newTrack.setTrackParameters(ts.getParameters(), B); newTrack.setCovarianceMatrix(new SymmetricMatrix(5, ts.getCovMatrix(), true)); + if(debug)System.out.println("Track state at perigee : " +ts.toString()); } // Add the hits to the track int firstHit_idx = -1; @@ -640,17 +676,71 @@ 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, saveTrackStateAtIntercept); if (ts != null) newTrack.getTrackStates().add(ts); + if(debug)System.out.println(ts.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); + /* 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 + 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; + + /* + + 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 = 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()); + */ + //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); + + 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()); + // 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); + + Vec targetFace=origin; + Plane targetPlane = new Plane(targetFace, new Vec(0., 1., 0.)); + HelixState helixAtTarget=kT.getHelixAtPlane(targetPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane + 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); + if(debug)System.out.println("Helix at Target from helix.toTrackState"); + if(debug)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); @@ -658,6 +748,7 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { newTrack.getTrackStates().add(ts_target); } } + */ // other track properties newTrack.setChisq(kT.chi2); @@ -759,6 +850,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); @@ -1345,7 +1437,7 @@ public void compareAllTracks(String trackCollectionName, EventHeader event, Arra TrackState ts1 = null; for (TrackState state : tkr.getTrackStates()) { //System.out.format("Track state %d: location=%d\n", tkr.getTrackStates().indexOf(state), state.getLocation()); - if (state.getLocation() == TrackState.AtIP) { + if (state.getLocation() == TrackState.AtPerigee) { ts1 = state; break; } @@ -1567,7 +1659,7 @@ public void plotGBLtracks(String path, EventHeader event) { printWriter3.format("$tkr%d << EOD\n", tracksGBL.indexOf(tkr)); for (TrackState state : tkr.getTrackStates()) { int loc = state.getLocation(); - if (loc != state.AtIP && loc != state.AtCalorimeter && loc != state.AtOther && loc != state.AtVertex) { + if (loc != state.AtPerigee && loc != state.AtCalorimeter && loc != state.AtOther && loc != state.AtVertex) { double [] pnt = state.getReferencePoint(); printWriter3.format(" %10.6f %10.6f %10.6f\n", pnt[0], pnt[2], -pnt[1]); } 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/KalmanPatRecPlots.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java index a8446d7740..bce8d2e4bc 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java @@ -851,7 +851,7 @@ void process(EventHeader event, List sensors, ArrayList[] } List stLst = tkrGBL.getTrackStates(); for (TrackState st : stLst) { - if (st.getLocation() == TrackState.AtIP) { + if (st.getLocation() == TrackState.AtPerigee) { double d0 = st.getParameter(0); aida.histogram1D("GBL d0").fill(d0); double z0 = st.getParameter(3); 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"); }