diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java b/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java index d7c48631e..277d8f4cb 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java @@ -48,13 +48,21 @@ public double getBeamEnergy() { protected boolean debug = false; protected boolean outputPlots = false; protected String outputPlotDir = "DQMOutputPlots/"; - + //define all of the collection strings here so the the variable names are all consistent + //makes it easier to do steering for multiple drivers. + protected String finalStateParticlesColName = "FinalStateParticles"; + protected String unconstrainedV0CandidatesColName = "UnconstrainedV0Candidates"; + protected String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates"; + protected String targetConV0CandidatesColName = "TargetConstrainedV0Candidates"; + @Override protected void detectorChanged(Detector detector) { BeamEnergyCollection beamEnergyCollection = this.getConditionsManager() .getCachedConditions(BeamEnergyCollection.class, "beam_energies").getCachedData(); - if (beamEnergy == null && beamEnergyCollection != null && beamEnergyCollection.size() != 0) + if (beamEnergy == null && beamEnergyCollection != null && beamEnergyCollection.size() != 0){ beamEnergy = beamEnergyCollection.get(0).getBeamEnergy(); + System.out.println(this.getClass().getName()+":: setting beamEnergy = "+beamEnergy); + } else { LOGGER.log(Level.WARNING, "warning: beam energy not found. Using a 6.6 GeV as the default energy"); beamEnergy = 6.6; @@ -63,9 +71,26 @@ protected void detectorChanged(Detector detector) { is2019Run = DatabaseConditionsManager.isPhys2019Run(this.getConditionsManager().getRun()); } + public void setFinalStateParticlesColName(String fspColName){ + finalStateParticlesColName=fspColName; + } + + public void setUnconstrainedV0CandidatesColName(String ucV0ColName){ + unconstrainedV0CandidatesColName=ucV0ColName; + } + + public void setBeamConV0CandidatesColName(String bscV0ColName){ + beamConV0CandidatesColName=bscV0ColName; + } + + public void setTargetConV0CandidatesColName(String tcV0ColName){ + targetConV0CandidatesColName=tcV0ColName; + } + String triggerType = "all";// allowed types are "" (blank) or "all", singles0, singles1, pairs0,pairs1 public boolean isGBL = false; - + public boolean isKF = false; + public void setTriggerType(String type) { this.triggerType = type; } @@ -73,6 +98,9 @@ public void setTriggerType(String type) { public void setIsGBL(boolean isgbl) { this.isGBL = isgbl; } + public void setIsKF(boolean iskf) { + this.isKF = iskf; + } public void setRecoVersion(String recoVersion) { this.recoVersion = recoVersion; diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/FinalStateMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/FinalStateMonitoring.java index 70f55b104..b8a350f0f 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/FinalStateMonitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/FinalStateMonitoring.java @@ -41,8 +41,6 @@ public class FinalStateMonitoring extends DataQualityMonitor { private static Logger LOGGER = Logger.getLogger(FinalStateMonitoring.class.getPackage().getName()); - String finalStateParticlesColName = "FinalStateParticles"; - String[] fpQuantNames = {"nEle_per_Event", "nPos_per_Event", "nPhoton_per_Event", "nUnAssociatedTracks_per_Event", "avg_delX_at_ECal", "avg_delY_at_ECal", "avg_E_Over_P", "avg_mom_beam_elec", "sig_mom_beam_elec"}; // some counters @@ -92,22 +90,20 @@ public class FinalStateMonitoring extends DataQualityMonitor { /* number of unassocaited tracks/event */ IHistogram1D nUnAssTracksHisto; - public void setFinalStateParticlesColName(String fsp) { - this.finalStateParticlesColName = fsp; - } - @Override protected void detectorChanged(Detector detector) { super.detectorChanged(detector); double maxFactor = 1.5; double feeMomentumCut = 0.75; // this number, multiplied by the beam energy, is the actual cut - beamEnergy = 4.5; - System.out.println("Using beamEnergy = " + beamEnergy); + System.out.println(this.getClass().getName()+":: setting beamEnergy = "+beamEnergy); LOGGER.info("Setting up the plotter"); aida.tree().cd("/"); String trkType = "SeedTrack/"; if (isGBL) trkType = "GBLTrack/"; + if (isKF) + trkType = "KFTrack/"; + /* Final State Particle Quantities */ /* plot electron & positron momentum separately */ @@ -247,12 +243,6 @@ public void process(EventHeader event) { if (fsCluster == null) throw new RuntimeException("isPhoton==true but no cluster found: should never happen"); double ene = fsPart.getEnergy(); - // TODO: mg-May 14, 2014....I would like to do this!!!! - // double xpos = fsCluster.getPositionAtShowerMax(false)[0];// false-->assume a photon instead of - // electron from calculating shower depth - // double ypos = fsCluster.getPositionAtShowerMax(false)[1]; - // but I can't because ReconParticles don't know about HPSEcalClusters, and casting it as one doesn't - // seem to work Hep3Vector clusterPosition = new BasicHep3Vector(fsCluster.getPosition()); double xpos = clusterPosition.x(); double ypos = clusterPosition.y(); @@ -271,9 +261,9 @@ public void process(EventHeader event) { double ene = fsPart.getEnergy(); double eOverP = ene / mom.magnitude(); Hep3Vector clusterPosition = new BasicHep3Vector(fsCluster.getPosition());// this gets position at - // shower max assuming it's an - // electron/positron - Hep3Vector trackPosAtEcal = TrackUtils.extrapolateTrack(fsTrack, clusterPosition.z()); + //get this from stored track state...don't do extrapolation by hand (commented out below) + Hep3Vector trackPosAtEcal=new BasicHep3Vector(TrackUtils.getTrackStateAtECal(fsTrack).getReferencePoint()); + // Hep3Vector trackPosAtEcal = TrackUtils.extrapolateTrack(fsTrack, clusterPosition.z()); double dx = trackPosAtEcal.x() - clusterPosition.x();// remember track vs detector coords double dy = trackPosAtEcal.y() - clusterPosition.y();// remember track vs detector coords diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/KFTrackingMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/KFTrackingMonitoring.java new file mode 100644 index 000000000..de8f30d68 --- /dev/null +++ b/analysis/src/main/java/org/hps/analysis/dataquality/KFTrackingMonitoring.java @@ -0,0 +1,581 @@ +package org.hps.analysis.dataquality; + +import hep.aida.IFitFactory; +import hep.aida.IFitResult; +import hep.aida.IFitter; +import hep.aida.IHistogram1D; +import hep.aida.IHistogram2D; +import hep.physics.vec.BasicHep3Matrix; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; +import java.util.ArrayList; +import java.util.Collection; +import java.util.List; +import java.util.Map; +import java.util.logging.Logger; +import org.hps.recon.tracking.CoordinateTransformations; +import org.hps.recon.tracking.TrackUtils; +import org.hps.recon.tracking.gbl.GBLKinkData; +import org.lcsim.detector.tracker.silicon.HpsSiSensor; +import org.lcsim.event.EventHeader; +import org.lcsim.event.GenericObject; +import org.lcsim.event.LCRelation; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.RelationalTable; +import org.lcsim.event.Track; +import org.lcsim.event.TrackerHit; +import org.lcsim.fit.helicaltrack.HelixUtils; +import org.lcsim.geometry.Detector; +import org.lcsim.util.aida.AIDA; + +/** + * DQM driver for reconstructed track quantities plots things like number of + * tracks/event, chi^2, track parameters (d0/z0/theta/phi/curvature) + */ +// TODO: Add some quantities for DQM monitoring: e.g. , , etc +public class KFTrackingMonitoring extends DataQualityMonitor { + + private static Logger LOGGER = Logger.getLogger(KFTrackingMonitoring.class.getPackage().getName()); + + private String trackCollectionName = "KalmanFullTracks"; + private final String trackerName = "Tracker"; + private static final String nameStrip = "Tracker_TestRunModule_"; + String ecalSubdetectorName = "Ecal"; + String ecalCollectionName = "EcalClustersCorr"; + private Detector detector = null; + private List sensors; + + int nEvents = 0; + int nTotTracks = 0; + int nTotHits = 0; + double sumd0 = 0; + double sumz0 = 0; + double sumslope = 0; + double sumchisq = 0; + private final String plotDir = "Tracks/"; + private final String positronDir = "Positrons/"; + private final String electronDir = "Electrons/"; + private final String topDir = "Top/"; + private final String botDir = "Bottom/"; + private final String hthplotDir = "TrackHits/"; + private final String timeresidDir = "HitTimeResiduals/"; + private final String kinkDir = "Kinks/"; + String[] trackingQuantNames = {"avg_N_tracks", "avg_N_hitsPerTrack", "avg_d0", "avg_z0", "avg_absslope", "avg_chi2"}; + int nmodules = 7; + IHistogram1D[] hthTop = new IHistogram1D[nmodules]; + IHistogram1D[] hthBot = new IHistogram1D[nmodules]; + IHistogram2D[] xvsyTop = new IHistogram2D[nmodules]; + IHistogram2D[] xvsyBot = new IHistogram2D[nmodules]; + IHistogram2D[] xvsyOnTrackTop = new IHistogram2D[nmodules]; + IHistogram2D[] xvsyOnTrackBot = new IHistogram2D[nmodules]; + IHistogram2D[] timevstimeTop = new IHistogram2D[nmodules]; + IHistogram2D[] timevstimeBot = new IHistogram2D[nmodules]; + IHistogram2D[] timevstimeOnTrackTop = new IHistogram2D[nmodules]; + IHistogram2D[] timevstimeOnTrackBot = new IHistogram2D[nmodules]; + IHistogram1D[] deltaTOnTrackTop = new IHistogram1D[nmodules]; + IHistogram1D[] deltaTOnTrackBot = new IHistogram1D[nmodules]; + + IHistogram1D trkChi2; + IHistogram1D nTracks; + IHistogram1D trkd0; + IHistogram1D trkphi; + IHistogram1D trkomega; + IHistogram1D trklam; + IHistogram1D trkz0; + IHistogram1D nHits; + IHistogram1D trackMeanTime; + IHistogram1D trackRMSTime; + IHistogram2D trackNhitsVsChi2; + IHistogram2D trackChi2RMSTime; + IHistogram1D seedRMSTime; + IHistogram2D trigTimeTrackTime; + IHistogram1D trigTime; + + IHistogram1D trkXAtECALTop; + IHistogram1D trkXAtECALBot; + IHistogram1D trkYAtECALTop; + IHistogram1D trkYAtECALBot; + IHistogram2D trkAtECAL; + + IHistogram1D trkChi2Pos; + IHistogram1D trkChi2Ele; + IHistogram1D trkChi2Top; + IHistogram1D trkChi2Bot; + + IHistogram1D nTracksPos; + IHistogram1D nTracksEle; + IHistogram1D nTracksTop; + IHistogram1D nTracksBot; + + IHistogram1D trkd0Pos; + IHistogram1D trkd0Ele; + IHistogram1D trkd0Top; + IHistogram1D trkd0Bot; + + IHistogram1D trkphiPos; + IHistogram1D trkphiEle; + IHistogram1D trkphiTop; + IHistogram1D trkphiBot; + + IHistogram1D trkomegaPos; + IHistogram1D trkomegaEle; + IHistogram1D trkomegaTop; + IHistogram1D trkomegaBot; + + IHistogram1D trklamPos; + IHistogram1D trklamEle; + IHistogram1D trklamTop; + IHistogram1D trklamBot; + + IHistogram1D trkz0Pos; + IHistogram1D trkz0Ele; + IHistogram1D trkz0Top; + IHistogram1D trkz0Bot; + + IHistogram1D nHitsPos; + IHistogram1D nHitsEle; + IHistogram1D nHitsTop; + IHistogram1D nHitsBot; + + IHistogram1D trkTimePos; + IHistogram1D trkTimeEle; + IHistogram1D trkTimeTop; + IHistogram1D trkTimeBot; + + IHistogram2D d0VsPhi0; + IHistogram2D d0Vsomega; + IHistogram2D d0Vslambda; + IHistogram2D d0Vsz0; + IHistogram2D phi0Vsomega; + IHistogram2D phi0Vslambda; + IHistogram2D phi0Vsz0; + IHistogram2D omegaVslambda; + IHistogram2D omegaVsz0; + IHistogram2D lamdbaVsz0; + + IHistogram2D chi2VsD0; + IHistogram2D chi2VsPhi0; + IHistogram2D chi2VsOmega; + IHistogram2D chi2VsLambda; + IHistogram2D chi2VsZ0; + + IHistogram2D beamAngleXY; + IHistogram2D beamAngleThetaPhi; + + IHistogram1D L1Iso; + IHistogram1D L12Iso; + IHistogram2D d0VsL1Iso; + IHistogram2D d0VsL12Iso; + + double d0Cut = 5.0; + double phiCut = 0.2; + double omegaCut = 0.0005; + double lambdaCut = 0.1; + double z0Cut = 1.0; + double timeCut = 24.0; + double maxChi2 = 100; + + double ecalXRange = 500; + double ecalYRange = 100; + + public void setTrackCollectionName(String trackCollectionName) { + this.trackCollectionName = trackCollectionName; + } + + @Override + protected void detectorChanged(Detector detector) { + super.detectorChanged(detector); + this.detector = detector; + aida.tree().cd("/"); + + trkChi2 = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track Chi2", 200, 0, maxChi2); + nTracks = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Tracks per Event", 6, 0, 6); + trkd0 = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0", 100, -d0Cut, d0Cut); + trkphi = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "sinphi", 100, -phiCut, phiCut); + trkomega = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "omega", 100, -omegaCut, omegaCut); + trklam = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "tan(lambda)", 100, -lambdaCut, lambdaCut); + trkz0 = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "z0", 100, -z0Cut, z0Cut); + nHits = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Hits per Track", 12, 3, 15); + trackMeanTime = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Mean time of hits on track", 100, -timeCut, timeCut); + trackRMSTime = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "RMS time of hits on track", 100, 0., 15.); + trackChi2RMSTime = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track chi2 vs. RMS time of hits", 100, 0., 15., 25, 0, maxChi2); + trackNhitsVsChi2 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track nhits vs. chi2", 100, 0, maxChi2, 10, 5, 15); + trigTimeTrackTime = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Trigger phase vs. mean time of hits", 100, -timeCut, timeCut, 6, 0, 24.0); + trigTime = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Trigger phase", 6, 0., 24.); + seedRMSTime = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "RMS time of hits on seed layers", 100, 0., 15.); + trkXAtECALTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track X at ECAL: Top", 100, -ecalXRange, ecalXRange); + trkXAtECALBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track X at ECAL: Bot", 100, -ecalXRange, ecalXRange); + trkYAtECALTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track Y at ECAL: Top", 100, 0, ecalYRange); + trkYAtECALBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track Y at ECAL: Bot", 100, 0, ecalYRange); + trkAtECAL = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track XY at ECAL", 100, -ecalXRange, ecalXRange, 100, -ecalYRange, ecalYRange); + for (int i = 0; i < nmodules; i++) { + double maxHTHX = 150.0; + double maxHTHY = 50.0; + if (i < 2) { + maxHTHX = 10; + maxHTHY = 15; + } else if (i < 4) { + maxHTHX = 50; + maxHTHY = 30; + } + xvsyTop[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Top", 250, -maxHTHX, maxHTHX, 30, 0, maxHTHY); + xvsyBot[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom", 250, -maxHTHX, maxHTHX, 30, 0, maxHTHY); + xvsyOnTrackTop[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Top, Hits On Tracks", 250, -maxHTHX, maxHTHX, 30, 0, maxHTHY); + xvsyOnTrackBot[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom, Hits On Tracks", 250, -maxHTHX, maxHTHX, 30, 0, maxHTHY); + timevstimeTop[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Top: Hit Times", 30, -15, 15, 30, -15, 15); + timevstimeBot[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom: Hit Times", 30, -15, 15, 30, -15, 15); + hthTop[i] = aida.histogram1D(hthplotDir + triggerType + "/" + "Module " + i + " Top: Track Hits", 25, 0, 25); + hthBot[i] = aida.histogram1D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom: Track Hits", 25, 0, 25); + timevstimeOnTrackTop[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Top: Hit Times, Hits On Tracks", 30, -15, 15, 30, -15, 15); + timevstimeOnTrackBot[i] = aida.histogram2D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom: Hit Times, Hits On Tracks", 30, -15, 15, 30, -15, 15); + deltaTOnTrackTop[i] = aida.histogram1D(hthplotDir + triggerType + "/" + "Module " + i + " Top: Hit Time Differences, Hits On Tracks", 50, -25, 25); + deltaTOnTrackBot[i] = aida.histogram1D(hthplotDir + triggerType + "/" + "Module " + i + " Bottom: Hit Time Differences, Hits On Tracks", 50, -25, 25); + } + + trkChi2Pos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "Track Chi2", 50, 0, maxChi2); + nTracksPos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "Tracks per Event", 6, 0, 6); + trkd0Pos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "d0", 100, -d0Cut, d0Cut); + trkphiPos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "sinphi", 100, -phiCut, phiCut); + trkomegaPos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "omega", 100, -omegaCut, omegaCut); + trklamPos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "tan(lambda)", 100, -lambdaCut, lambdaCut); + trkz0Pos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "z0", 100, -z0Cut, z0Cut); + nHitsPos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "Hits per Track", 12, 3, 15); + trkTimePos = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + positronDir + "Mean time of hits on track", 100, -timeCut, timeCut); + + trkChi2Ele = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "Track Chi2", 50, 0, maxChi2); + nTracksEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "Tracks per Event", 6, 0, 6); + trkd0Ele = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "d0", 100, -d0Cut, d0Cut); + trkphiEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "sinphi", 100, -phiCut, phiCut); + trkomegaEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "omega", 100, -omegaCut, omegaCut); + trklamEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "tan(lambda)", 100, -lambdaCut, lambdaCut); + trkz0Ele = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "z0", 100, -z0Cut, z0Cut); + nHitsEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "Hits per Track", 12, 3, 15); + trkTimeEle = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + electronDir + "Mean time of hits on track", 100, -timeCut, timeCut); + + trkChi2Top = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "Track Chi2", 50, 0, maxChi2); + nTracksTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "Tracks per Event", 6, 0, 6); + trkd0Top = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "d0", 100, -d0Cut, d0Cut); + trkphiTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "sinphi", 100, -phiCut, phiCut); + trkomegaTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "omega", 100, -omegaCut, omegaCut); + trklamTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "tan(lambda)", 100, 0.0, lambdaCut); + trkz0Top = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "z0", 100, -z0Cut, z0Cut); + nHitsTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "Hits per Track", 12, 3, 15); + trkTimeTop = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + topDir + "Mean time of hits on track", 100, -timeCut, timeCut); + + trkChi2Bot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "Track Chi2",50, 0, maxChi2); + nTracksBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "Tracks per Event", 6, 0, 6); + trkd0Bot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "d0", 100, -d0Cut, d0Cut); + trkphiBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "sinphi", 100, -phiCut, phiCut); + trkomegaBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "omega", 100, -omegaCut, omegaCut); + trklamBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "tan(lambda)", 100, 0, lambdaCut); + trkz0Bot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "z0", 100, -z0Cut, z0Cut); + nHitsBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "Hits per Track", 12, 3, 15); + trkTimeBot = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + botDir + "Mean time of hits on track", 100, -timeCut, timeCut); + + //correlation plots + d0VsPhi0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs sinphi", 50, -d0Cut, d0Cut, 50, -phiCut, phiCut); + d0Vsomega = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs omega", 50, -d0Cut, d0Cut, 50, -omegaCut, omegaCut); + d0Vslambda = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs tan(lambda)", 50, -d0Cut, d0Cut, 50, -lambdaCut, lambdaCut); + d0Vsz0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs z0", 50, -d0Cut, d0Cut, 50, -z0Cut, z0Cut); + + phi0Vsomega = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "phi0 vs omega", 50, -phiCut, phiCut, 50, -omegaCut, omegaCut); + phi0Vslambda = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "phi0 vs tan(lambda)", 50, -phiCut, phiCut, 50, -lambdaCut, lambdaCut); + phi0Vsz0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "phi0 vs z0", 50, -phiCut, phiCut, 50, -z0Cut, z0Cut); + + omegaVslambda = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "omega vs tan(lambda)", 50, -omegaCut, omegaCut, 50, -lambdaCut, lambdaCut); + omegaVsz0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "omega vs z0", 50, -omegaCut, omegaCut, 50, -z0Cut, z0Cut); + + lamdbaVsz0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "lambda vs z0", 50, -lambdaCut, lambdaCut, 50, -z0Cut, z0Cut); + + chi2VsD0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "chi2 vs d0", 50, -d0Cut, d0Cut, 50, 0.0, maxChi2); + chi2VsPhi0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "chi2 vs sinphi", 50, -phiCut, phiCut, 50, 0.0, maxChi2); + chi2VsOmega = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "chi2 vs omega", 50, -omegaCut, omegaCut, 50, 0.0, maxChi2); + chi2VsLambda = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "chi2 vs lambda", 50, -lambdaCut, lambdaCut, 50, 0.0, maxChi2); + chi2VsZ0 = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "chi2 vs z0", 50, -z0Cut, z0Cut, 50, 0.0, 50.0); + + beamAngleXY = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "angles around beam axis: theta_y vs theta_x", 100, -0.1, 0.1, 100, -0.1, 0.1); + beamAngleThetaPhi = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "angles around beam axis: theta vs phi", 100, -Math.PI, Math.PI, 100, 0, 0.25); + + L1Iso = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "L1 isolation", 100, -5.0, 5.0); + L12Iso = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "L1-2 isolation", 100, -5.0, 5.0); + d0VsL1Iso = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs L1 isolation", 100, -5.0, 5.0, 50, -d0Cut, d0Cut); + d0VsL12Iso = aida.histogram2D(plotDir + trackCollectionName + "/" + triggerType + "/" + "d0 vs L1-2 isolation", 100, -5.0, 5.0, 50, -d0Cut, d0Cut); + + // Make a list of SiSensors in the SVT. + sensors = this.detector.getSubdetector(trackerName).getDetectorElement().findDescendants(HpsSiSensor.class); + + // Setup the occupancy plots. + aida.tree().cd("/"); + for (HpsSiSensor sensor : sensors) { + //IHistogram1D occupancyPlot = aida.histogram1D(sensor.getName().replaceAll("Tracker_TestRunModule_", ""), 640, 0, 639); + IHistogram1D hitTimeResidual = PlotAndFitUtilities.createSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + timeresidDir + "hitTimeResidual_", sensor, 100, -20, 20,true); + IHistogram1D lambdaKink = PlotAndFitUtilities.createSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "lambdaKink_", sensor, 100, -5e-3, 5e-3,true); + IHistogram1D phiKink = PlotAndFitUtilities.createSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "phiKink_", sensor, 100, -5e-3, 5e-3,true); + IHistogram2D lambdaKink2D = PlotAndFitUtilities.createSensorPlot2D(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "lambdaKinkVsOmega_", sensor, 100, -omegaCut, omegaCut, 100, -5e-3, 5e-3,true); + IHistogram2D phiKink2D = PlotAndFitUtilities.createSensorPlot2D(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "phiKinkVsOmega_", sensor, 100, -omegaCut, omegaCut, 100, -5e-3, 5e-3,true); + } + } + + @Override + public void process(EventHeader event) { + + int runNumber = event.getRunNumber(); + + aida.tree().cd("/"); + + + + //check to see if this event is from the correct trigger (or "all"); + if (!matchTrigger(event)) + return; + + + if (!event.hasCollection(Track.class, trackCollectionName)) { + nTracks.fill(0); + return; + } + nEvents++; + List tracks = event.get(Track.class, trackCollectionName); + nTotTracks += tracks.size(); + nTracks.fill(tracks.size()); + int cntEle = 0; + int cntPos = 0; + int cntTop = 0; + int cntBot = 0; + for (Track trk : tracks) { + Hep3Vector trackPosAtEcalFace = TrackUtils.getTrackPositionAtEcal(trk, runNumber); + double xAtECal = trackPosAtEcalFace.x(); + double yAtECal = trackPosAtEcalFace.y(); + if (yAtECal > 0) { + trkXAtECALTop.fill(xAtECal); + trkYAtECALTop.fill(yAtECal); + } else { + trkXAtECALBot.fill(xAtECal); + trkYAtECALBot.fill(Math.abs(yAtECal)); + } + trkAtECAL.fill(xAtECal, yAtECal); + nTotHits += trk.getTrackerHits().size(); + + double d0 = trk.getTrackStates().get(0).getD0(); + double sinphi0 = Math.sin(trk.getTrackStates().get(0).getPhi()); + double omega = trk.getTrackStates().get(0).getOmega(); + double lambda = trk.getTrackStates().get(0).getTanLambda(); + double z0 = trk.getTrackStates().get(0).getZ0(); + trkChi2.fill(trk.getChi2()); + nHits.fill(trk.getTrackerHits().size()); + trackNhitsVsChi2.fill(trk.getChi2(), trk.getTrackerHits().size()); + trkd0.fill(trk.getTrackStates().get(0).getD0()); + trkphi.fill(Math.sin(trk.getTrackStates().get(0).getPhi())); + trkomega.fill(trk.getTrackStates().get(0).getOmega()); + trklam.fill(trk.getTrackStates().get(0).getTanLambda()); + trkz0.fill(trk.getTrackStates().get(0).getZ0()); + d0VsPhi0.fill(d0, sinphi0); + d0Vsomega.fill(d0, omega); + d0Vslambda.fill(d0, lambda); + d0Vsz0.fill(d0, z0); + phi0Vsomega.fill(sinphi0, omega); + phi0Vslambda.fill(sinphi0, lambda); + phi0Vsz0.fill(sinphi0, z0); + omegaVslambda.fill(omega, lambda); + omegaVsz0.fill(omega, z0); + lamdbaVsz0.fill(lambda, z0); + chi2VsD0.fill(d0, trk.getChi2()); + chi2VsPhi0.fill(sinphi0, trk.getChi2()); + chi2VsOmega.fill(omega, trk.getChi2()); + chi2VsLambda.fill(lambda, trk.getChi2()); + chi2VsZ0.fill(z0, trk.getChi2()); + + BasicHep3Matrix rot = new BasicHep3Matrix(); + rot.setActiveEuler(Math.PI / 2, -0.0305, -Math.PI / 2); + Hep3Vector dir = CoordinateTransformations.transformVectorToDetector(HelixUtils.Direction(TrackUtils.getHTF(trk), 0)); + Hep3Vector dirRotated = VecOp.mult(rot, dir); + + double beamPhi = Math.atan2(dirRotated.y(), dirRotated.x()); + double beamTheta = Math.acos(dirRotated.z()); + + beamAngleXY.fill(dirRotated.x(), dirRotated.y()); + beamAngleThetaPhi.fill(beamPhi, beamTheta); + /* + Double[] isolations = TrackUtils.getIsolations(trk, hitToStrips, hitToRotated); + double l1Iso = Double.MAX_VALUE; + double l12Iso = Double.MAX_VALUE; + + for (int i = 0; i < isolations.length; i++) + if (isolations[i] != null) { + if (i == 0 || i == 1) + if (Math.abs(isolations[i]) < Math.abs(l1Iso)) + l1Iso = isolations[i]; + if (i < 4) + if (Math.abs(isolations[i]) < Math.abs(l12Iso)) + l12Iso = isolations[i]; + } + L1Iso.fill(l1Iso); + L12Iso.fill(l12Iso); + d0VsL1Iso.fill(l1Iso, d0); + d0VsL12Iso.fill(l12Iso, d0); + */ + int nStrips = 0; + double meanTime = 0; + double meanSeedTime = 0; + + List stripHits = new ArrayList(); + + for (TrackerHit hit : trk.getTrackerHits()) { + + stripHits.add(hit); + int stripLayer = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getLayerNumber(); + nStrips++; + meanTime += hit.getTime(); + + } + meanTime /= nStrips; + + double rmsTime = 0; + double rmsSeedTime = 0; + + stripHits = TrackUtils.sortHits(stripHits); + for (TrackerHit hts : stripHits) { + rmsTime += Math.pow(hts.getTime() - meanTime, 2); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) hts.getRawHits().get(0)).getDetectorElement(); + int layer = sensor.getLayerNumber(); + if (layer <= 6) + rmsSeedTime += Math.pow(hts.getTime() - meanSeedTime, 2); + PlotAndFitUtilities.getSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + timeresidDir + "hitTimeResidual_", sensor,true).fill((hts.getTime() - meanTime) * nStrips / (nStrips - 1)); //correct residual for bias + } + rmsTime = Math.sqrt(rmsTime / nStrips); + trackMeanTime.fill(meanTime); + trackRMSTime.fill(rmsTime); + trackChi2RMSTime.fill(rmsTime, trk.getChi2()); + trigTimeTrackTime.fill(meanTime, event.getTimeStamp() % 24); + trigTime.fill(event.getTimeStamp() % 24); + + GenericObject kinkData = GBLKinkData.getKinkData(event, trk); + if (kinkData != null) + for (int i = 0; i < stripHits.size(); i++) { + TrackerHit hts = stripHits.get(i); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) hts.getRawHits().get(0)).getDetectorElement(); +// int layer = sensor.getLayerNumber(); + double lambdaKink = GBLKinkData.getLambdaKink(kinkData, i); + double phiKink = GBLKinkData.getPhiKink(kinkData, i); +// System.out.format("%d %d %f %f\n", i, layer, lambdaKink, phiKink); + + PlotAndFitUtilities.getSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "lambdaKink_", sensor,true).fill(lambdaKink); + PlotAndFitUtilities.getSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "phiKink_", sensor,true).fill(phiKink); + PlotAndFitUtilities.getSensorPlot2D(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "lambdaKinkVsOmega_", sensor,true).fill(trk.getTrackStates().get(0).getOmega(), lambdaKink); + PlotAndFitUtilities.getSensorPlot2D(plotDir + trackCollectionName + "/" + triggerType + "/" + kinkDir + "phiKinkVsOmega_", sensor,true).fill(trk.getTrackStates().get(0).getOmega(), phiKink); + } + + if (trk.getTrackStates().get(0).getOmega() < 0) {//positrons + trkChi2Pos.fill(trk.getChi2()); + nHitsPos.fill(trk.getTrackerHits().size()); + trkd0Pos.fill(trk.getTrackStates().get(0).getD0()); + trkphiPos.fill(Math.sin(trk.getTrackStates().get(0).getPhi())); + trkomegaPos.fill(trk.getTrackStates().get(0).getOmega()); + trklamPos.fill(trk.getTrackStates().get(0).getTanLambda()); + trkz0Pos.fill(trk.getTrackStates().get(0).getZ0()); + trkTimePos.fill(meanTime); + cntPos++; + } else { + trkChi2Ele.fill(trk.getChi2()); + nHitsEle.fill(trk.getTrackerHits().size()); + trkd0Ele.fill(trk.getTrackStates().get(0).getD0()); + trkphiEle.fill(Math.sin(trk.getTrackStates().get(0).getPhi())); + trkomegaEle.fill(trk.getTrackStates().get(0).getOmega()); + trklamEle.fill(trk.getTrackStates().get(0).getTanLambda()); + trkz0Ele.fill(trk.getTrackStates().get(0).getZ0()); + trkTimeEle.fill(meanTime); + cntEle++; + } + + if (trk.getTrackStates().get(0).getTanLambda() < 0) { + trkChi2Bot.fill(trk.getChi2()); + nHitsBot.fill(trk.getTrackerHits().size()); + trkd0Bot.fill(trk.getTrackStates().get(0).getD0()); + trkphiBot.fill(Math.sin(trk.getTrackStates().get(0).getPhi())); + trkomegaBot.fill(trk.getTrackStates().get(0).getOmega()); + trklamBot.fill(Math.abs(trk.getTrackStates().get(0).getTanLambda())); + trkz0Bot.fill(trk.getTrackStates().get(0).getZ0()); + trkTimeBot.fill(meanTime); + + cntBot++; + } else { + trkChi2Top.fill(trk.getChi2()); + nHitsTop.fill(trk.getTrackerHits().size()); + trkd0Top.fill(trk.getTrackStates().get(0).getD0()); + trkphiTop.fill(Math.sin(trk.getTrackStates().get(0).getPhi())); + trkomegaTop.fill(trk.getTrackStates().get(0).getOmega()); + trklamTop.fill(Math.abs(trk.getTrackStates().get(0).getTanLambda())); + trkz0Top.fill(trk.getTrackStates().get(0).getZ0()); + trkTimeTop.fill(meanTime); + cntTop++; + } + + sumd0 += trk.getTrackStates().get(0).getD0(); + sumz0 += trk.getTrackStates().get(0).getZ0(); + sumslope += Math.abs(trk.getTrackStates().get(0).getTanLambda()); + sumchisq += trk.getChi2(); + +// System.out.format("%d seed strips, RMS time %f\n", nSeedStrips, rmsSeedTime); +// System.out.format("%d strips, mean time %f, RMS time %f\n", nStrips, meanTime, rmsTime); + } + nTracksTop.fill(cntTop); + nTracksBot.fill(cntBot); + nTracksPos.fill(cntPos); + nTracksEle.fill(cntEle); + } + + @Override + public void calculateEndOfRunQuantities() { + IFitFactory fitFactory = AIDA.defaultInstance().analysisFactory().createFitFactory(); + IFitter fitter = fitFactory.createFitter("chi2"); + + for (HpsSiSensor sensor : sensors) { + //IHistogram1D occupancyPlot = aida.histogram1D(sensor.getName().replaceAll("Tracker_TestRunModule_", ""), 640, 0, 639); + IHistogram1D hitTimeResidual = PlotAndFitUtilities.getSensorPlot(plotDir + trackCollectionName + "/" + triggerType + "/" + timeresidDir + "hitTimeResidual_", sensor,true); + IFitResult result = fitGaussian(hitTimeResidual, fitter, "range=\"(-20.0,20.0)\""); + if (result != null) + System.out.format("%s\t%f\t%f\t%d\t%d\t%f\n", getNiceSensorName(sensor), result.fittedParameters()[1], result.fittedParameters()[2], sensor.getFebID(), sensor.getFebHybridID(), sensor.getT0Shift()); + } + + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[0], (double) nTotTracks / nEvents); + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[1], (double) nTotHits / nTotTracks); + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[2], sumd0 / nTotTracks); + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[3], sumz0 / nTotTracks); + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[4], sumslope / nTotTracks); + monitoredQuantityMap.put(trackCollectionName + " " + triggerType + " " + trackingQuantNames[5], sumchisq / nTotTracks); + } + + IFitResult fitGaussian(IHistogram1D h1d, IFitter fitter, String range) { + double[] init = {h1d.maxBinHeight(), h1d.mean(), h1d.rms()}; + IFitResult ifr = null; + try { + ifr = fitter.fit(h1d, "g", init, range); + } catch (RuntimeException ex) { + LOGGER.info(this.getClass().getSimpleName() + ": caught exception in fitGaussian"); + } + return ifr; +// double[] init = {20.0, 0.0, 1.0, 20, -1}; +// return fitter.fit(h1d, "g+p1", init, range); + } + + @Override + public void printDQMData() { + LOGGER.info("ReconMonitoring::printDQMData"); + for (Map.Entry entry : monitoredQuantityMap.entrySet()) + LOGGER.info(entry.getKey() + " = " + entry.getValue()); + } + + @Override + public void printDQMStrings() { + for (Map.Entry entry : monitoredQuantityMap.entrySet()) + LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;"); + } + + private String getNiceSensorName(HpsSiSensor sensor) { + return sensor.getName().replaceAll(nameStrip, "") + .replace("module", "mod") + .replace("layer", "lyr") + .replace("sensor", "sens"); + } + +} diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/MollerMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/MollerMonitoring.java index 144b31b08..65ca02e87 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/MollerMonitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/MollerMonitoring.java @@ -75,11 +75,6 @@ public String toString() { // "Isolation"}; // private int firstVertexingCut = 0; - private final String finalStateParticlesColName = "FinalStateParticles"; - private final String unconstrainedV0CandidatesColName = "UnconstrainedMollerCandidates"; - // private final String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates"; - // private final String targetV0ConCandidatesColName = "TargetConstrainedV0Candidates"; - // private final String trackListName = "MatchedTracks"; private final String[] fpQuantNames = {"nV0_per_Event", "avg_BSCon_mass", "avg_BSCon_Vx", "avg_BSCon_Vy", "avg_BSCon_Vz", "sig_BSCon_Vx", "sig_BSCon_Vy", "sig_BSCon_Vz", "avg_BSCon_Chi2"}; @@ -355,13 +350,6 @@ public void process(EventHeader event) { if (!event.hasCollection(ReconstructedParticle.class, unconstrainedV0CandidatesColName)) { return; } - // if (!event.hasCollection(ReconstructedParticle.class, beamConV0CandidatesColName)) { - // return; - // } - // if (!event.hasCollection(ReconstructedParticle.class, targetV0ConCandidatesColName)) { - // return; - // } - // check to see if this event is from the correct trigger (or "all"); if (!matchTrigger(event)) { return; diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/MuonCandidateMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/MuonCandidateMonitoring.java index 46101cbb5..0aa4b273a 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/MuonCandidateMonitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/MuonCandidateMonitoring.java @@ -27,10 +27,6 @@ public class MuonCandidateMonitoring extends DataQualityMonitor { private static Logger LOGGER = Logger.getLogger(V0Monitoring.class.getPackage().getName()); - private String finalStateParticlesColName = "FinalStateParticles"; - private String unconstrainedV0CandidatesColName = "UnconstrainedV0Candidates"; - private String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates"; - private String targetV0ConCandidatesColName = "TargetConstrainedV0Candidates"; private String clusterCollectionName = "EcalClustersCorr"; private int nRecoEvents = 0; diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java index 8481d8d8e..044755f9c 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java @@ -61,7 +61,7 @@ public class TrackingMonitoring extends DataQualityMonitor { private final String electronDir = "Electrons/"; private final String topDir = "Top/"; private final String botDir = "Bottom/"; - private final String hthplotDir = "HelicalTrackHits/"; + private final String hthplotDir = "TrackHits/"; private final String timeresidDir = "HitTimeResiduals/"; private final String kinkDir = "Kinks/"; String[] trackingQuantNames = {"avg_N_tracks", "avg_N_hitsPerTrack", "avg_d0", "avg_z0", "avg_absslope", "avg_chi2"}; @@ -194,6 +194,7 @@ public void setTrackCollectionName(String trackCollectionName) { protected void detectorChanged(Detector detector) { super.detectorChanged(detector); this.detector = detector; + aida.tree().cd("/"); trkChi2 = aida.histogram1D(plotDir + trackCollectionName + "/" + triggerType + "/" + "Track Chi2", 200, 0, maxChi2); diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java index 935a98ad8..e05c556c2 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/TridentMonitoring.java @@ -86,11 +86,7 @@ public String toString() { private static final int TRIDENT = 0; private static final int VERTEX = 1; - private final String finalStateParticlesColName = "FinalStateParticles"; - private final String unconstrainedV0CandidatesColName = "UnconstrainedV0Candidates"; - // private final String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates"; - // private final String targetV0ConCandidatesColName = "TargetConstrainedV0Candidates"; - // private final String trackListName = "MatchedTracks"; + private final String[] fpQuantNames = {"nV0_per_Event", "avg_BSCon_mass", "avg_BSCon_Vx", "avg_BSCon_Vy", "avg_BSCon_Vz", "sig_BSCon_Vx", "sig_BSCon_Vy", "sig_BSCon_Vz", "avg_BSCon_Chi2"}; @@ -219,14 +215,15 @@ public String toString() { private final IHistogram2D[][] cutVertexZVsMass = new IHistogram2D[Cut.nCuts][2]; private final double plotsMinMass = 0.01; - private final double plotsMaxMass = 0.1; + private final double plotsMaxMass = 0.3; // clean up event first - private final int nTrkMax = 5; - private final int nPosMax = 1; + private final int nTrkMax = 10; + private final int nPosMax = 5; private final double maxChi2SeedTrack = 7.0; private double maxChi2GBLTrack = 100.0; + private double maxChi2KFTrack = 100.0; private double maxUnconVertChi2 = 100.0; private double maxBsconVertChi2 = 1000.0; // disable by default @@ -251,12 +248,12 @@ public String toString() { private double v0BsconVxCut = 20.0; // disable by default // track quality cuts - private final double beamPCut = 1.0; + private final double beamPCut = 1.0;//times the beam energy private final double minPCut = 0.1; // private double trkPyMax = 0.2; // private double trkPxMax = 0.2; private final double radCut = 0.6; - private final double trkTimeDiff = 10.0; + private final double trkTimeDiff = 20.0; private final double clusterTimeDiffCut = 5.0; private double l1IsoMin = 0.250; @@ -328,7 +325,7 @@ public void setBeamPosY(double beamPosY) { protected void detectorChanged(Detector detector) { super.detectorChanged(detector); /* tab */ - LOGGER.info("TridendMonitoring::detectorChanged Setting up the plotter"); + LOGGER.info("TridentMonitoring::detectorChanged Setting up the plotter"); beamAxisRotation.setActiveEuler(Math.PI / 2, -0.0305, -Math.PI / 2); beamEnergy = 4.5; aida.tree().cd("/"); @@ -689,7 +686,7 @@ public void process(EventHeader event) { stdDev += Math.pow(time - mean, 2); stdDev /= (hitTimes.size() - 1); stdDev = Math.sqrt(stdDev); - + /* Double[] eleIso = TrackUtils.getIsolations(eleTrack, hitToStrips, hitToRotated); Double[] posIso = TrackUtils.getIsolations(posTrack, hitToStrips, hitToRotated); double minPositiveIso = 9999; @@ -734,7 +731,7 @@ public void process(EventHeader event) { } double minIso = Math.min(Math.abs(minPositiveIso), Math.abs(minNegativeIso)); - + */ double tEle = TrackUtils.getTrackTime(eleTrack, hitToStrips, hitToRotated); double tPos = TrackUtils.getTrackTime(posTrack, hitToStrips, hitToRotated); @@ -770,8 +767,10 @@ public void process(EventHeader event) { // start applying cuts EnumSet bits = EnumSet.noneOf(Cut.class); // System.out.println("Starting to apply cuts"); +// boolean trackQualityCut = Math.max(electron.getTracks().get(0).getChi2(), positron.getTracks().get(0) +// .getChi2()) < (isGBL ? maxChi2GBLTrack : maxChi2SeedTrack); boolean trackQualityCut = Math.max(electron.getTracks().get(0).getChi2(), positron.getTracks().get(0) - .getChi2()) < (isGBL ? maxChi2GBLTrack : maxChi2SeedTrack); + .getChi2()) < maxChi2KFTrack; // System.out.println("Track quality cuts: " + trackQualityCut); if (trackQualityCut) bits.add(Cut.TRK_QUALITY); @@ -837,11 +836,13 @@ public void process(EventHeader event) { // System.out.println("frontHitsCut: "+frontHitsCut); if (frontHitsCut) bits.add(Cut.FRONT_HITS); - boolean isoCut = minIso > l1IsoMin; + // boolean isoCut = minIso > l1IsoMin; // System.out.println("isoCut: "+isoCut); - if (!frontHitsCut || isoCut) // diagnostic plots look better if failing the front hits cut makes you pass +// if (!frontHitsCut || isoCut) // diagnostic plots look better if failing the front hits cut makes you pass // // this one - bits.add(Cut.ISOLATION); + boolean isoCut=true; + if(isoCut) + bits.add(Cut.ISOLATION); for (Cut cut : Cut.values()) if (bits.contains(cut)) { // if (cut.ordinal() == Cut.firstVertexingCut)// if we get here, we've passed all non-vertexing cuts @@ -858,10 +859,10 @@ public void process(EventHeader event) { // && uncV0.getMomentum().magnitude() > radCut * beamEnergy) if (1 == 1) switch (cut) { - case ISOLATION: - l1Iso.fill(minIso); - zVsL1Iso.fill(minIso, v0Vtx.z()); - break; + // case ISOLATION: + // l1Iso.fill(minIso); + // zVsL1Iso.fill(minIso, v0Vtx.z()); + // break; case EVENT_QUALITY: eventTrkCount.fill(ntrk); eventPosCount.fill(npos); @@ -1072,14 +1073,14 @@ public void process(EventHeader event) { @Override // TODO: Change from System.out to use logger instead. public void printDQMData() { - System.out.println("TridendMonitoring::printDQMData"); + System.out.println("TridentMonitoring::printDQMData"); for (Entry entry : monitoredQuantityMap.entrySet()) System.out.println(entry.getKey() + " = " + entry.getValue()); System.out.println("*******************************"); - System.out.println("TridendMonitoring::Tridend Selection Summary: " + (isGBL ? "GBLTrack" : "SeedTrack")); + System.out.println("TridentMonitoring::Trident Selection Summary: " + (isGBL ? "GBLTrack" : (isKF ? "KalmanTrack" : "SeedTrack"))); - System.out.println("\t\t\tTridend Selection Summary"); + System.out.println("\t\t\tTrident Selection Summary"); System.out .println("******************************************************************************************"); System.out.println(String.format("Number of V0:\t%8.0f\t%8.6f\t%8.6f\t%8.6f\n", nRecoV0, diff --git a/analysis/src/main/java/org/hps/analysis/dataquality/V0Monitoring.java b/analysis/src/main/java/org/hps/analysis/dataquality/V0Monitoring.java index 724435fc4..13ed7d7d5 100644 --- a/analysis/src/main/java/org/hps/analysis/dataquality/V0Monitoring.java +++ b/analysis/src/main/java/org/hps/analysis/dataquality/V0Monitoring.java @@ -41,10 +41,7 @@ public class V0Monitoring extends DataQualityMonitor { private static Logger LOGGER = Logger.getLogger(V0Monitoring.class.getPackage().getName()); - private String finalStateParticlesColName = "FinalStateParticles"; - private String unconstrainedV0CandidatesColName = "UnconstrainedV0Candidates"; - private String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates"; - private String targetV0ConCandidatesColName = "TargetConstrainedV0Candidates"; + private String[] fpQuantNames = {"nV0_per_Event", "avg_BSCon_mass", "avg_BSCon_Vx", "avg_BSCon_Vy", "avg_BSCon_Vz", "sig_BSCon_Vx", "sig_BSCon_Vy", "sig_BSCon_Vz", "avg_BSCon_Chi2"}; // some counters @@ -162,12 +159,7 @@ public class V0Monitoring extends DataQualityMonitor { @Override protected void detectorChanged(Detector detector) { super.detectorChanged(detector); - - BeamEnergyCollection beamEnergyCollection = this.getConditionsManager().getCachedConditions(BeamEnergyCollection.class, "beam_energies").getCachedData(); - beamEnergy = beamEnergyCollection.get(0).getBeamEnergy(); - - System.out.println("Using beamEnergy = " + beamEnergy); - + feeMomentumCut = 0.75 * beamEnergy; // GeV v0ESumMinCut = 0.8 * beamEnergy; @@ -178,6 +170,8 @@ protected void detectorChanged(Detector detector) { molPSumMax = 1.25 * beamEnergy; beambeamCut = 0.80 * beamEnergy; + double maxChi2=100.0; + beamAxisRotation.setActiveEuler(Math.PI / 2, -0.0305, -Math.PI / 2); LOGGER.info("Setting up the plotter"); @@ -187,6 +181,9 @@ protected void detectorChanged(Detector detector) { if (isGBL) trkType = "GBLTrack/"; + if (isKF) + trkType = "KFTrack/"; + double maxMass = .1 * beamEnergy; double maxMassMoller = .1 * Math.sqrt(beamEnergy); @@ -202,13 +199,13 @@ protected void detectorChanged(Detector detector) { unconVz = aida.histogram1D(plotDir + trkType + triggerType + "/" + unconstrainedV0CandidatesColName + "/" + "Vz (mm)", 50, -50, 50); unconESum = aida.histogram1D(plotDir + trkType + triggerType + "/" + unconstrainedV0CandidatesColName + "/" - + "ESum (GeV)", 50, 0.1 * beamEnergy, v0ESumMaxCut); + + "PSum (GeV)", 50, 0.1 * beamEnergy, v0ESumMaxCut); unconChi2 = aida.histogram1D(plotDir + trkType + triggerType + "/" + unconstrainedV0CandidatesColName + "/" + "Chi2", 25, 0, 25); unconVzVsChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + unconstrainedV0CandidatesColName + "/" + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50); unconChi2VsTrkChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + unconstrainedV0CandidatesColName - + "/" + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25); + + "/" + "Chi2 vs. total track chi2", 50, 0, maxChi2, 50, 0, maxChi2); /* beamspot constrained */ bsconMass = aida.histogram1D(plotDir + trkType + triggerType + "/" + beamConV0CandidatesColName + "/" + "Mass (GeV)", 100, 0, maxMass); @@ -221,28 +218,28 @@ protected void detectorChanged(Detector detector) { bsconESum = aida.histogram1D( plotDir + trkType + triggerType + "/" + beamConV0CandidatesColName + "/" + "ESum (GeV)", 50, 0.1 * beamEnergy, v0ESumMaxCut); bsconChi2 = aida.histogram1D(plotDir + trkType + triggerType + "/" + beamConV0CandidatesColName + "/" + "Chi2", - 25, 0, 25); + 50, 0, maxChi2); bsconVzVsChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + beamConV0CandidatesColName + "/" - + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50); + + "Vz vs. Chi2", 50, 0, maxChi2, 50, -50, 50); bsconChi2VsTrkChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + beamConV0CandidatesColName + "/" - + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25); + + "Chi2 vs. total track chi2", 50, 0, maxChi2, 50, 0, maxChi2); /* target constrained */ - tarconMass = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" + tarconMass = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + "Mass (GeV)", 100, 0, maxMass); - tarconVx = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" + tarconVx = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + "Vx (mm)", 50, -1, 1); - tarconVy = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" + tarconVy = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + "Vy (mm)", 50, -1, 1); - tarconVz = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" + tarconVz = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + "Vz (mm)", 50, -10, 10); - tarconESum = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" + tarconESum = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + "ESum (GeV)", 50, 0.1 * beamEnergy, v0ESumMaxCut); - tarconChi2 = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" - + "Chi2", 25, 0, 25); - tarconVzVsChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName + "/" - + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50); - tarconChi2VsTrkChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + targetV0ConCandidatesColName - + "/" + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25); + tarconChi2 = aida.histogram1D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + + "Chi2", 50, 0, maxChi2); + tarconVzVsChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + "/" + + "Vz vs. Chi2", 50, 0, maxChi2, 50, -50, 50); + tarconChi2VsTrkChi2 = aida.histogram2D(plotDir + trkType + triggerType + "/" + targetConV0CandidatesColName + + "/" + "Chi2 vs. total track chi2", 50, 0, maxChi2, 50, 0, maxChi2); nV0 = aida .histogram1D(plotDir + trkType + triggerType + "/" + xtra + "/" + "Number of V0 per event", 10, 0, 10); @@ -380,8 +377,8 @@ public void process(EventHeader event) { throw new IllegalArgumentException("Collection missing: " + unconstrainedV0CandidatesColName); if (!event.hasCollection(ReconstructedParticle.class, beamConV0CandidatesColName)) throw new IllegalArgumentException("Collection missing: " + beamConV0CandidatesColName); - if (!event.hasCollection(ReconstructedParticle.class, targetV0ConCandidatesColName)) - throw new IllegalArgumentException("Collection missing: " + targetV0ConCandidatesColName); + if (!event.hasCollection(ReconstructedParticle.class, targetConV0CandidatesColName)) + throw new IllegalArgumentException("Collection missing: " + targetConV0CandidatesColName); // check to see if this event is from the correct trigger (or "all"); if (!matchTrigger(event)) @@ -405,7 +402,7 @@ public void process(EventHeader event) { unconVx.fill(vtxPosRot.x()); unconVy.fill(vtxPosRot.y()); unconVz.fill(vtxPosRot.z()); - unconESum.fill(uncV0.getEnergy()); + unconESum.fill(uncV0.getMomentum().magnitude()); unconMass.fill(uncV0.getMass()); unconChi2.fill(uncVert.getChi2()); unconVzVsChi2.fill(uncVert.getChi2(), vtxPosRot.z()); @@ -453,6 +450,7 @@ public void process(EventHeader event) { Math.max(uncV0.getParticles().get(0).getTracks().get(0).getChi2(), uncV0.getParticles().get(1) .getTracks().get(0).getChi2()), uncVert.getPosition().z()); + /* Double[] eleIso = TrackUtils.getIsolations(ele.getTracks().get(0), hitToStrips, hitToRotated); Double[] posIso = TrackUtils.getIsolations(pos.getTracks().get(0), hitToStrips, hitToRotated); if (eleIso[0] != null && posIso[0] != null) { @@ -461,7 +459,7 @@ public void process(EventHeader event) { double minL1Iso = Math.min(eleL1Iso, posL1Iso); VtxZVsL1Iso.fill(minL1Iso, uncVert.getPosition().z()); } - + */ double pe = ele.getMomentum().magnitude(); double pp = pos.getMomentum().magnitude(); Hep3Vector pEleRot = VecOp.mult(beamAxisRotation, ele.getMomentum()); @@ -504,7 +502,7 @@ public void process(EventHeader event) { bsconVx.fill(vtxPosRot.x()); bsconVy.fill(vtxPosRot.y()); bsconVz.fill(vtxPosRot.z()); - bsconESum.fill(bsV0.getEnergy()); + bsconESum.fill(bsV0.getMomentum().magnitude()); bsconMass.fill(bsV0.getMass()); bsconChi2.fill(bsVert.getChi2()); bsconVzVsChi2.fill(bsVert.getChi2(), vtxPosRot.z()); @@ -519,7 +517,7 @@ public void process(EventHeader event) { } nV0.fill(v0Count); List targetConstrainedV0List = event.get(ReconstructedParticle.class, - targetV0ConCandidatesColName); + targetConV0CandidatesColName); for (ReconstructedParticle tarV0 : targetConstrainedV0List) { if (isGBL != TrackType.isGBL(tarV0.getType())) @@ -530,7 +528,7 @@ public void process(EventHeader event) { tarconVx.fill(vtxPosRot.x()); tarconVy.fill(vtxPosRot.y()); tarconVz.fill(vtxPosRot.z()); - tarconESum.fill(tarV0.getEnergy()); + tarconESum.fill(tarV0.getMomentum().magnitude()); tarconMass.fill(tarV0.getMass()); tarconChi2.fill(tarVert.getChi2()); tarconVzVsChi2.fill(tarVert.getChi2(), vtxPosRot.z()); @@ -597,6 +595,7 @@ else if (!hasSharedStrips(ele1, fsPart, hitToStrips, hitToRotated)) } // look at "Moller" events (if that's what they really are + /* if (ele1.getMomentum().magnitude() + ele2.getMomentum().magnitude() > molPSumMin && ele1.getMomentum().magnitude() + ele2.getMomentum().magnitude() < molPSumMax && (p1.magnitude() < beambeamCut && p2.magnitude() < beambeamCut)) { @@ -647,6 +646,7 @@ else if (!hasSharedStrips(ele1, fsPart, hitToStrips, hitToRotated)) pEleVsthetaMoller.fill(p2.magnitude(), theta2); thetaEleVsthetaMoller.fill(theta1, theta2); } + */ } } diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed.lcsim index eca765dd9..49849d4da 100644 --- a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed.lcsim +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed.lcsim @@ -175,7 +175,7 @@ false false v0skim_parameters_ver1.txt - ${outputFile}.slcio + ${outputFile}_v0skim.slcio diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed_dataqual.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed_dataqual.lcsim new file mode 100644 index 000000000..cb74bcad4 --- /dev/null +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass2_recon_skimmed_dataqual.lcsim @@ -0,0 +1,236 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1000 + + + svt_readout_overlap_good + + + + + + + + CONFIG + + + + + + + + + WARNING + EcalClusters + + + EcalClusters + EcalClustersCorr + + + + SVTRawTrackerHits + + + .5 + 1 + Pileup + Migrad + false + true + false + false + true + true + false + true + false + + + 24.0 + 3.0 + false + 400 + 4.0 + 1.0 + 3.0 + 3.0 + true + true + false + + + 3 + 1 + 8.757651 + 38.0487 + 3.98915 + 11.777395 + 0 + 3 + 3 + 39.95028 + 8.186345 + + 13.52662 + 7.00678 + + 9.771546584 + 1.7652935 + 5 + false + 466 + .725912 + + 0.0 + 0.02 + 0.0 + 0.05 + 0.0 + 1.0 + 7.204329 + false + + + EcalClustersCorr + KalmanFullTracks + KalmanFullTracks + TrackClusterMatcherMinDistance + UnconstrainedV0Candidates_KF + UnconstrainedV0Vertices_KF + BeamspotConstrainedV0Candidates_KF + BeamspotConstrainedV0Vertices_KF + TargetConstrainedV0Candidates_KF + TargetConstrainedV0Vertices_KF + FinalStateParticles_KF + true + false + false + 0.180 + 0.350 + 0.100 + 0.100 + 0.900 + 7.0 + 7.0 + 0.0 + 40.0 + 40 + 40 + false + true + true + false + true + UnconstrainedMollerCandidates_KF + UnconstrainedMollerVertices_KF + BeamspotConstrainedMollerCandidates_KF + BeamspotConstrainedMollerVertices_KF + TargetConstrainedMollerCandidates_KF + TargetConstrainedMollerVertices_KF + + + true + false + false + false + v0skim_parameters_ver1.txt + ${outputFile}_v0skim.slcio + + + + + all + EcalClustersCorr + false + + + all + + + all + KalmanFullTracks + + + FinalStateParticles_KF + all + true + + + FinalStateParticles_KF + UnconstrainedV0Candidates_KF + BeamspotConstrainedV0Candidates_KF + TargetConstrainedV0Candidates_KF + all + true + + + FinalStateParticles_KF + UnconstrainedV0Candidates_KF + BeamspotConstrainedV0Candidates_KF + TargetConstrainedV0Candidates_KF + all + true + + + ${outputFile}_data_quality_plots.root + + + + + FPGAData HelicalTrackHitRelations HelicalTrackHits HelicalTrackMCRelations KFGBLStripClusterData KFGBLStripClusterDataRelations ReadoutTimestamps RotatedHelicalTrackHitRelations RotatedHelicalTrackHits RotatedHelicalTrackMCRelations SVTFittedRawTrackerHits SVTShapeFitParameters SVTTrueHitRelations StripClusterer_SiTrackerHitStrip1D SVTRawTrackerHits FADCGenericHits HodoReadoutHits HodoCalHits EcalReadoutHits EcalUncalHits HodoGenericClusters VTPBank EcalClusters + ${outputFile}.slcio + + + + 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 4537c2d00..71ef5d4b0 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java +++ b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java @@ -1539,18 +1539,26 @@ public static RelationalTable getHitToRotatedTable(EventHeader event) { public static double getTrackTime(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) { double meanTime = 0; - List stripHits = getStripHits(track, hitToStrips, hitToRotated); + List stripHits=new ArrayList(); + // if(hitToStrips == null || hitToRotated == null) + // stripHits= track.getTrackerHits(); + //else + stripHits = getStripHits(track, hitToStrips, hitToRotated); + for (TrackerHit hit : stripHits) { meanTime += hit.getTime(); } meanTime /= stripHits.size(); return meanTime; - } + } public static double getTrackTimeSD(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) { double meanTime = getTrackTime(track, hitToStrips, hitToRotated); - List stripHits = getStripHits(track, hitToStrips, hitToRotated); - + List stripHits=new ArrayList(); + // if(hitToStrips == null || hitToRotated == null) + // stripHits= track.getTrackerHits(); + //else + stripHits = getStripHits(track, hitToStrips, hitToRotated); double sdTime = 0; for (TrackerHit hit : stripHits) { sdTime += Math.pow(meanTime - hit.getTime(), 2); @@ -1561,14 +1569,18 @@ public static double getTrackTimeSD(Track track, RelationalTable hitToStrips, Re } public static List getStripHits(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) { - List hits = new ArrayList(); - for (TrackerHit hit : track.getTrackerHits()) { - hits.addAll(hitToStrips.allFrom(hitToRotated.from(hit))); - } + List hits = new ArrayList(); + if(hitToStrips == null || hitToRotated == null) + hits=track.getTrackerHits(); + else + for (TrackerHit hit : track.getTrackerHits()) + hits.addAll(hitToStrips.allFrom(hitToRotated.from(hit))); + return hits; } + public static List sortHits(Collection hits) { List hitList = new ArrayList(hits); Collections.sort(hitList, new LayerComparator()); @@ -1593,14 +1605,20 @@ public int compare(TrackerHit o1, TrackerHit o2) { * @return */ public static int numberOfSharedStrips(Track track1, Track track2, RelationalTable hitToStrips, RelationalTable hitToRotated) { - Set track1hits = new HashSet(getStripHits(track1, hitToStrips, hitToRotated)); - int nShared = 0; - for (TrackerHit hit : track2.getTrackerHits()) { - for (TrackerHit hts : (Set) hitToStrips.allFrom(hitToRotated.from(hit))) { - if (track1hits.contains(hts)) { - nShared++; - } - } + Set track1hits =null; + Set track2hits =null; + // if(hitToStrips == null || hitToRotated == null){ //this is a kalman track + // track1hits = new HashSet(track1.getTrackerHits()); + // track2hits = new HashSet(track2.getTrackerHits()); + //} else { + track1hits=new HashSet(getStripHits(track1, hitToStrips, hitToRotated)); + track2hits=new HashSet(getStripHits(track2, hitToStrips, hitToRotated)); + //} + int nShared = 0; + for (TrackerHit hit : track2hits) { + if (track1hits.contains(hit)) { + nShared++; + } } return nShared; } @@ -1622,7 +1640,7 @@ public static boolean hasSharedStrips(Track track1, Track track2, RelationalTabl return true; } } - + public static int getLayer(TrackerHit strip) { if (strip == null) { System.out.println("Strip is null?????");