From 69b74affdb458ea9c39379c706cfe04790d01e5f Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Sat, 5 Jul 2025 07:43:56 -0700 Subject: [PATCH] add tests for kalman tracking; FEE, Moller, and V0 2016 --- .../it/PhysRun2016FeeReconKalmanTest.java | 216 ++++++++++++ .../it/PhysRun2016MollerReconKalmanTest.java | 312 ++++++++++++++++++ .../test/it/PhysRun2016V0ReconKalmanTest.java | 300 +++++++++++++++++ 3 files changed, 828 insertions(+) create mode 100644 integration-tests/src/test/java/org/hps/test/it/PhysRun2016FeeReconKalmanTest.java create mode 100644 integration-tests/src/test/java/org/hps/test/it/PhysRun2016MollerReconKalmanTest.java create mode 100644 integration-tests/src/test/java/org/hps/test/it/PhysRun2016V0ReconKalmanTest.java 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 + } + } +}