From 7585594e214f7f768b73a20f0b84c245673d7b30 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 20 Sep 2024 09:35:43 -0700 Subject: [PATCH 1/4] make mctweaking directory and move files there. Make StripHitSmearer.java copied from the killer. This may not even build yet --- .../KalmanSlopeBasedTrackHitKiller.java | 0 .../SlopeBasedTrackHitKiller.java | 0 .../{ => mctweaking}/StripHitKiller.java | 0 .../tracking/mctweaking/StripHitSmearer.java | 545 ++++++++++++++++++ 4 files changed, 545 insertions(+) rename tracking/src/main/java/org/hps/recon/tracking/{ => mctweaking}/KalmanSlopeBasedTrackHitKiller.java (100%) rename tracking/src/main/java/org/hps/recon/tracking/{ => mctweaking}/SlopeBasedTrackHitKiller.java (100%) rename tracking/src/main/java/org/hps/recon/tracking/{ => mctweaking}/StripHitKiller.java (100%) create mode 100644 tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java diff --git a/tracking/src/main/java/org/hps/recon/tracking/KalmanSlopeBasedTrackHitKiller.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/KalmanSlopeBasedTrackHitKiller.java similarity index 100% rename from tracking/src/main/java/org/hps/recon/tracking/KalmanSlopeBasedTrackHitKiller.java rename to tracking/src/main/java/org/hps/recon/tracking/mctweaking/KalmanSlopeBasedTrackHitKiller.java diff --git a/tracking/src/main/java/org/hps/recon/tracking/SlopeBasedTrackHitKiller.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/SlopeBasedTrackHitKiller.java similarity index 100% rename from tracking/src/main/java/org/hps/recon/tracking/SlopeBasedTrackHitKiller.java rename to tracking/src/main/java/org/hps/recon/tracking/mctweaking/SlopeBasedTrackHitKiller.java diff --git a/tracking/src/main/java/org/hps/recon/tracking/StripHitKiller.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitKiller.java similarity index 100% rename from tracking/src/main/java/org/hps/recon/tracking/StripHitKiller.java rename to tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitKiller.java diff --git a/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java new file mode 100644 index 000000000..7bffec4a3 --- /dev/null +++ b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java @@ -0,0 +1,545 @@ +package org.hps.recon.tracking; + +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStream; +import java.io.InputStreamReader; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.logging.Level; +import java.util.logging.Logger; +import java.util.regex.Matcher; +import java.util.regex.Pattern; +import org.lcsim.detector.tracker.silicon.ChargeCarrier; +import org.lcsim.detector.tracker.silicon.HpsSiSensor; +import org.lcsim.detector.tracker.silicon.SiSensorElectrodes; +import org.lcsim.event.EventHeader; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.TrackerHit; +import org.lcsim.geometry.Detector; +import org.lcsim.lcio.LCIOUtil; +import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D; +import org.lcsim.util.Driver; + +/** + * + * @author mgraham created 5/14/19 + * This driver will remove 1d strip clusters from the + * "StripClusterer_SiTrackerHitStrip1D" (default) + * collection based on a channel-based (data/mc) ratio file. + * + * Careful..the names of the ratio files are important! Format + * should be: + * _LX_top/bottom_stereo/axial_slot/hole.txt + * + * mg...this only works for L1 and L6 smearing at the moment + * ...unfortunately some of these things are hard coded + */ +public class StripHitSmearer extends Driver { + + //IMPORTANT...the layer, top/bottom/stereo/axial/slot/hole are derived from these names!!! + Set ratioFiles = new HashSet(); + private List _sensorsToSmear = new ArrayList(); + //List of Sensors + private List sensors = null; + private static final String SUBDETECTOR_NAME = "Tracker"; + private static Pattern layerPattern = Pattern.compile("L(\\d+)(t|b)"); + private String stripHitInputCollectionName = "StripClusterer_SiTrackerHitStrip1D"; + private boolean _debug = false; + // instead of just removing hit at a given channel, remove all hits within Nsigma + private boolean removeHitsWithinNSig = false; + // used sigma-weighted ratios + private boolean useSigmaWeightedRatios = false; + double nSig = 5; + double sigmaUL1 = 0.1; //100 microns...this is roughly the 5-hit track projection error in U for layer 1, plotted in SVTHitLevelPlots + double sigmaUL6 = 0.25; //250 microns...this is roughly the 5-hit track projection error in U for layer 6*1.5 because pull is too wide, plotted in SVTHitLevelPlots + // double sigmaUL6 = 0.05; //50 microns...this is narrow just to check the effect. + // double firstSensorSmearFactor=3.0; + /// double secondSensorSmearFactor=2.0; + double firstSensorSmearFactor = 1.0; + double secondSensorSmearFactor = 1.0; + //these are just used for debugging... + int checkHitsChannel = 4; + int checkHitsTotal = 0; + int checkHitsPassed = 0; + double checkHitsRatio = 0; + + double maxLayer = 12; + /// + + private List siClusters = new ArrayList(); + + private Map _siClustersAcceptMap = new HashMap(); + private Map _finalSiClustersAcceptMap = new HashMap(); + + private Map _smearingFractionsL1 = new HashMap(); + private Map _smearingFractionsL6 = new HashMap(); + + public void setRatioFiles(String[] ratioNames) { + System.out.println("Setting ratio files!!! " + ratioNames[0]); + this.ratioFiles = new HashSet(Arrays.asList(ratioNames)); + } + + public void setDebug(boolean debug) { + this._debug = debug; + } + + public void setRemoveHitsWithinNSig(boolean removeHitsWithinNSig) { + this.removeHitsWithinNSig = removeHitsWithinNSig; + } + + public void setUseSigmaWeightedRatios(boolean useSigmaWeightedRatios) { + this.useSigmaWeightedRatios = useSigmaWeightedRatios; + } + + public void setL1Sigma(double l1sigma) { + this.sigmaUL1 = l1sigma; + } + + public void setL6Sigma(double l6sigma) { + this.sigmaUL1 = l6sigma; + } + + public void setNSigma(double nSig) { + this.nSig = nSig; + } + + public StripHitSmearer() { + } + + @Override + public void detectorChanged(Detector detector) { + // Get the HpsSiSensor objects from the tracker detector element + sensors = detector.getSubdetector(SUBDETECTOR_NAME) + .getDetectorElement().findDescendants(HpsSiSensor.class); + // If the detector element had no sensors associated with it, throw + // an exception + if (sensors.size() == 0) + throw new RuntimeException("No sensors were found in this detector."); + + //parse the ratio names and register sensors to kill + String delims = "_|\\.";// this will split strings between one "_" or "." + for (String ratioFile : ratioFiles) { + System.out.println("StripHitSmearer::Using this ratioFile: " + ratioFile); + int layer = -1; + boolean top = false; + boolean stereo = false; + boolean slot = false; + System.out.println("Parsing ratioFile = " + ratioFile); + String[] tokens = ratioFile.split(delims); + Matcher m = layerPattern.matcher(tokens[1]); + if (m.find()) { + layer = Integer.parseInt(m.group(1)); + if (m.group(2).matches("t")) + top = true; + } else { + System.out.println("Couldn't find layer number!!! " + ratioFile); + continue; + } +// if (tokens[2].matches("top")) +// top = true; +// System.out.println(tokens[2]+" "+tokens[3]+" "+tokens[4]); + if (tokens[2].matches("stereo")) + stereo = true; + if (tokens[3].matches("slot")) + slot = true; + System.out.println("StripHitSmearer::Smearing this: " + + "layer = " + layer + "; top = " + top + "; stereo = " + stereo + + "; slot = " + slot); + this.registerSensor(layer, top, stereo, slot, ratioFile); + } + System.out.println("filling smearing fractions"); + _smearingFractionsL1 = this.fillSmearedFractions(sigmaUL1, 1); + _smearingFractionsL6 = this.fillSmearedFractions(sigmaUL6, 12); + } + + @Override + public void process(EventHeader event) { + // System.out.println("In process of SVTHitSmearer"); + + int nhitsremoved = 0; + _siClustersAcceptMap.clear(); + if (event.hasItem(stripHitInputCollectionName)) + siClusters = (List) event.get(stripHitInputCollectionName); + else { + System.out.println("StripHitSmearer::process No Input Collection Found?? " + stripHitInputCollectionName); + return; + } + if (_debug) + System.out.println("Number of SiClusters Found = " + siClusters.size()); + int oldClusterListSize = siClusters.size(); + for (TrackerHit siCluster : siClusters) { + boolean passHit = true; + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) siCluster.getRawHits().get(0)).getDetectorElement(); + for (SensorToSmear sensorToSmear : _sensorsToSmear) + if (sensorToSmear.matchSensor(sensor)) { + //ok, get hit channel and kill or not + Hep3Vector pos = globalToSensor(toHep3(siCluster.getPosition()), sensor); + int chan = getChan(pos, sensor); + + double ratio = 0; + if (useSigmaWeightedRatios) + ratio = this.getSmearedRatio(chan, sensorToSmear); + else + ratio = sensorToSmear.getRatio(chan); + double killFactor = (1-ratio)*sensorToSmear.getScaleSmearFactor(); + if (_debug) + System.out.println("Found a hit on a sensor to kill!!! " + sensor.getName() + " Layer = " + sensor.getLayerNumber() + + " isTop? " + sensorToSmear.getIsTop() + " isStereo? " + sensorToSmear.getIsStereo() + + " isSlot? " + sensorToSmear.getIsSlot() + " channel = " + chan + " ratio = " + ratio); + + ratio = 1-killFactor; + if (ratio != -666) { + double random = Math.random(); //throw a random number to see if this hit should be rejected + if (random > ratio) { + passHit = false; + if (_debug) + System.out.println("Smearing this hit Layer = " + sensor.getLayerNumber() + + " isTop? " + sensorToSmear.getIsTop() + "; isStereo? " + sensorToSmear.getIsStereo() + + " isSlot? " + sensorToSmear.getIsSlot() +" scaleSmearFactor= " + sensorToSmear.getScaleSmearFactor() + + " channel = " + chan + " kill ratio = " + ratio + " random = " + random); + nhitsremoved++; + } + } + if (chan == checkHitsChannel) { + checkHitsTotal++; + if (passHit) + checkHitsPassed++; + checkHitsRatio = ratio; + } + + } + _siClustersAcceptMap.put(siCluster, passHit); + } + +// if (_debug) + List tmpClusterList = getFinalHits(_siClustersAcceptMap); + // if (_debug) + if (_debug)System.out.println("New Cluster List Has " + tmpClusterList.size() + "; old List had " + oldClusterListSize); + if (_debug)System.out.println("number of hits removed for this event = " + nhitsremoved); + // TODO flag not used + int flag = LCIOUtil.bitSet(0, 31, true); // Turn on 64-bit cell ID. + event.remove(this.stripHitInputCollectionName); + event.put(this.stripHitInputCollectionName, tmpClusterList, SiTrackerHitStrip1D.class, 0, toString()); + + } + + @Override + public void endOfData() { + System.out.println("StripHitSmearer::endOfData channel = " + checkHitsChannel + + " ratio = " + checkHitsRatio + "; hits pass = " + checkHitsPassed + + " hits total = " + checkHitsTotal + " passRatio = " + ((double) checkHitsPassed) / checkHitsTotal); + } + + public void registerSensor(int layer, boolean isTop, boolean isStereo, boolean isSlot, String ratioFile) { + if (_debug) { + System.out.println("================ Registering New SensorToSmear =============== "); + System.out.println(" " + + "layer = " + layer + "; top = " + isTop + "; stereo = " + isStereo + + "; slot = " + isSlot + "; ratio file = " + ratioFile); + } + SensorToSmear newSensor = new SensorToSmear(layer, isTop, isStereo, isSlot, ratioFile); + System.out.println("newSensor isTop " + newSensor.getIsTop()); + if (newSensor.getIsTop() && !(newSensor.getIsStereo()) && newSensor.getLayer()==1) + newSensor.setScaleSmearFactor(firstSensorSmearFactor) ; + if (newSensor.getIsTop() && newSensor.getIsStereo() && newSensor.getLayer()==1) + newSensor.setScaleSmearFactor(secondSensorSmearFactor); + + if (!(newSensor.getIsTop()) && newSensor.getIsStereo() && newSensor.getLayer()==1) + newSensor.setScaleSmearFactor(firstSensorSmearFactor); + if (!(newSensor.getIsTop()) && !(newSensor.getIsStereo()) && newSensor.getLayer()==1) + newSensor.setScaleSmearFactor(secondSensorSmearFactor); + if (_debug) { + System.out.println(" channel 10 ratio for newSensor = "+newSensor.getRatio(10)); + System.out.println("================ Done With New SensorToSmear =============== "); + } + _sensorsToSmear.add(newSensor); + } + + //Return the HpsSiSensor for a given top/bottom track, layer, axial/stereo, and slot/hole + private HpsSiSensor getSensor(int layer, boolean isTop, boolean isAxial, boolean isHole) { + for (HpsSiSensor sensor : sensors) { + int senselayer = (sensor.getLayerNumber() + 1) / 2; + if (senselayer != layer) + continue; + if ((isTop && !sensor.isTopLayer()) || (!isTop && sensor.isTopLayer())) + continue; + if ((isAxial && !sensor.isAxial()) || (!isAxial && sensor.isAxial())) + continue; + if (layer < 4 && layer > 0) {//only for pre-2019 runs...fix this... + if (_debug)System.out.println("Matching with sensor = " + sensor.getName()); + return sensor; + } + else { + if ((!sensor.getSide().matches("ELECTRON") && isHole) || (sensor.getSide().matches("ELECTRON") && !isHole)) + continue; + + if (_debug)System.out.println("Matching with sensor = " + sensor.getName()); + return sensor; + } + } + return null; + } + + private List getFinalHits(Map _initialSiClustersAcceptMap) { + _finalSiClustersAcceptMap.clear(); + List tmpClusterList = new ArrayList(); + for (TrackerHit hit : _initialSiClustersAcceptMap.keySet()) { + boolean isAccept = _initialSiClustersAcceptMap.get(hit); + if (removeHitsWithinNSig && !isAccept) { // this hit got thrown out; lets throw out all hits within N-sigma + if (_debug) + System.out.println("Found a killed hit...now check for hits around it"); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement(); + Hep3Vector pos = globalToSensor(toHep3(hit.getPosition()), sensor); + int chan = getChan(pos, sensor); + double sigmaU = sigmaUL1; + if (sensor.getLayerNumber() > 10) //remember, layer number is sensor number (1-12) + sigmaU = sigmaUL6; + if (_debug) + System.out.println("Layer number = " + sensor.getLayerNumber() + " sigma = " + sigmaU); + Hep3Vector minPos = new BasicHep3Vector(pos.x() - nSig * sigmaU, pos.y(), pos.z()); + Hep3Vector maxPos = new BasicHep3Vector(pos.x() + nSig * sigmaU, pos.y(), pos.z()); + int minChan = getChan(minPos, sensor); + int maxChan = getChan(maxPos, sensor); + if (minChan > maxChan) { + int tmp = minChan; + minChan = maxChan; + maxChan = tmp; + } + if (_debug) + System.out.println("hit channel = " + chan + "; minChan = " + minChan + "; maxChan = " + maxChan); + //loop through initial clusters again + for (TrackerHit testhit : _initialSiClustersAcceptMap.keySet()) { + HpsSiSensor testsensor = (HpsSiSensor) ((RawTrackerHit) testhit.getRawHits().get(0)).getDetectorElement(); + if (sensor != testsensor) + continue; + Hep3Vector testpos = globalToSensor(toHep3(testhit.getPosition()), sensor); + int testchan = getChan(testpos, sensor); + if (_debug) + System.out.println("test channel = " + testchan); + if (testchan > maxChan || testchan < minChan) + continue; + //in range, set accept to false + if (_debug) + System.out.println("Found a hit within range of killed hit...killing it too!"); + _finalSiClustersAcceptMap.put(testhit, false); + } + _finalSiClustersAcceptMap.put(hit, isAccept); //add the original cluster to map (shouldn't matter since isAccept==false + } else if (!_finalSiClustersAcceptMap.containsKey(hit)) + _finalSiClustersAcceptMap.put(hit, isAccept); + } + //fill tmpClusterList with accepted hits + for (TrackerHit hit : _finalSiClustersAcceptMap.keySet()) + if (_finalSiClustersAcceptMap.get(hit)) + tmpClusterList.add(hit); + return tmpClusterList; + } + + private double getSmearedRatio(int seedChan, SensorToSmear sensorToSmear) { + double smearedRatio = 0; + double integralSum = 0; + Map _smearingFractions = new HashMap(); + if (sensorToSmear.getLayer() == 1) + _smearingFractions = _smearingFractionsL1; + else if (sensorToSmear.getLayer() == 6) + _smearingFractions = _smearingFractionsL6; + for (Map.Entry entry : _smearingFractions.entrySet()) { + int chan = seedChan + entry.getKey(); + if (chan < 0 || chan > sensorToSmear._sensor.getNumberOfChannels() - 1) + continue; + double chanIntegral = entry.getValue(); + double ratio = sensorToSmear.getRatio(chan); + if (ratio > 1.0 || ratio == -666) + ratio = 1.0; +// System.out.println("channel = " + chan + " ratio = " + // + ratio + "; smearedRatio = " + ratio * chanIntegral); + smearedRatio += ratio * chanIntegral; + integralSum += chanIntegral; + } + if (integralSum == 0) + smearedRatio = 1.0; + else + smearedRatio /= integralSum; + if (_debug) + System.out.println("getSmearedRatio: seed ratio = " + + sensorToSmear.getRatio(seedChan) + "; smearedRatio = " + smearedRatio); + return smearedRatio; + } + + private Map fillSmearedFractions(double sigmaU, int layer) { + Map _smearingFractions = new HashMap(); + double pitch = 0; + for (HpsSiSensor s : sensors) + if (s.getLayerNumber() == layer) + pitch = s.getReadoutStripPitch(); + if (pitch == 0) { + System.out.println("Couldn't find layer = " + layer); + return _smearingFractions; + } + + double normPitch = sigmaU / pitch; + int nStrips = (int) Math.floor(this.nSig * normPitch); + double totalInt = 0; + for (int i = -nStrips; i < 1; i++) {//just do negative bins (and 0)..positive side is same + double minIntLimit = (i - 0.5) / normPitch; + double maxIntLimit = (i + 0.5) / normPitch; + double minInt = this.computeGaussInt(minIntLimit, 1000); + double maxInt = this.computeGaussInt(maxIntLimit, 1000); + totalInt += maxInt - minInt; + _smearingFractions.put(i, maxInt - minInt); + System.out.println("_smearingFractions chan = " + i + " fraction = " + _smearingFractions.get(i)); + + } + for (int i = 1; i < nStrips + 1; i++) {//positive side is same + totalInt += _smearingFractions.get(-1 * i); + _smearingFractions.put(i, _smearingFractions.get(-1 * i)); + System.out.println("_smearingFractions chan = " + i + " fraction = " + _smearingFractions.get(i)); + + } + System.out.println("Layer = " + layer + "; totalInt = " + totalInt); + return _smearingFractions; + } + + public class SensorToSmear { + + int _layer = 1; + boolean _isStereo = false; + boolean _isSlot = false; + boolean _isTop = false; + String _ratioFile = "foobarTopL1Stereo.txt"; + HpsSiSensor _sensor = null; + double _scaleSmearFactor = 1.0; + Map _channelToRatioMap = new HashMap(); + + public SensorToSmear(int layer, boolean isTop, boolean isStereo, boolean isSlot, String ratioFile) { + if (_debug)System.out.println("Making new SensorToSmear layer = " + layer); + _layer = layer; + _isTop = isTop; + _isStereo = isStereo; + _isSlot = isSlot; + _ratioFile = ratioFile; + _sensor = getSensor(layer, isTop, !isStereo, !isSlot); + readRatioFile(); + } + + void setScaleSmearFactor(double scale) { + this._scaleSmearFactor=scale; + } + + int getLayer() { + return _layer; + } + + boolean getIsStereo() { + return _isStereo; + } + + boolean getIsSlot() { + return _isSlot; + } + + boolean getIsTop() { + return _isTop; + } + + String getRatioFile() { + return _ratioFile; + } + + double getScaleSmearFactor() { + return(this._scaleSmearFactor); + } + + boolean matchSensor(HpsSiSensor sensor) { + return _sensor == sensor; + } + + private void readRatioFile() { + String infile = "/org/hps/recon/tracking/efficiencyCorrections/" + _ratioFile; + InputStream inRatios = this.getClass().getResourceAsStream(infile); + if (_debug)System.out.println("StripHitSmearer::Reading ratio file " + infile); + BufferedReader reader = new BufferedReader(new InputStreamReader(inRatios)); + String line; + String delims = "[ ]+";// this will split strings between one or more spaces + try { + while ((line = reader.readLine()) != null) { + String[] tokens = line.split(delims); + if (_debug) System.out.println("channel number = " + tokens[0] + "; ratio = " + tokens[1]); + _channelToRatioMap.put(Integer.parseInt(tokens[0]), Double.parseDouble(tokens[1])); + } + } catch (IOException ex) { + Logger.getLogger(StripHitSmearer.class.getName()).log(Level.SEVERE, null, ex); + } + } + + public double getRatio(int channel) { + if (_debug)System.out.println("SensorToSmear:: getRatio: layer = " + this.getLayer() + "; channel " + channel); + if (_channelToRatioMap.get(channel) == null) + return -666; + return _channelToRatioMap.get(channel); + } + + } + + //Converts double array into Hep3Vector + private Hep3Vector toHep3(double[] arr) { + return new BasicHep3Vector(arr[0], arr[1], arr[2]); + } + //Returns channel number of a given position in the sensor frame + /* + private int getChan(Hep3Vector pos, HpsSiSensor sensor) { + double readoutPitch = sensor.getReadoutStripPitch(); + int nChan = sensor.getNumberOfChannels(); + double height = readoutPitch * nChan; + return (int) ((height / 2 - pos.x()) / readoutPitch); + } + */ + + private int getChan(Hep3Vector pos, HpsSiSensor sensor) { + SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + if (maxLayer>12 && sensor.getLayerNumber()<5) {//thin sensors + int row = electrodes.getRowNumber(pos); + int col = electrodes.getColumnNumber(pos); + return electrodes.getCellID(row,col); + }else{ + return electrodes.getCellID(pos); + } + } + //Converts position into sensor frame + private Hep3Vector globalToSensor(Hep3Vector trkpos, HpsSiSensor sensor) { + SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + if (electrodes == null) { + electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON); + System.out.println("Charge Carrier is NULL"); + } + return electrodes.getGlobalToLocal().transformed(trkpos); + } + + //Computes gaussian integral numerically from -inf to nSig + // this should be replaced by a simple error function but I can't find one + // in normal java libraries!!! + private double computeGaussInt(double nSig, int nSteps) { + double mean = 0; + double sigma = 1; + double dx = sigma * nSig / (double) nSteps; + double integral = 0; + for (int i = 0; i < nSteps; i++) { + double x = dx * (i + 0.5) + mean; + integral += dx * Gauss(x, mean, sigma); + } + return integral + 0.5; + } + + //Gaussian function + private double Gauss(double x, double mean, double sigma) { + return 1 / (Math.sqrt(2 * Math.PI * Math.pow(sigma, 2))) * Math.exp(-Math.pow(x - mean, 2) / (2 * Math.pow(sigma, 2))); + } + +} From 108a63ef374df896ea347d25117e19253ad902ad Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 19 Dec 2024 09:54:36 -0800 Subject: [PATCH 2/4] adding in time and position smearings plus an example steering file; this also includes some changes to KalmanPatRecDriver to related unbiased residuals to MP IDs instead of layer but i will separate this out to its own pull request as well at KFOutput --- ...Run2016FullReconMC_KF_Killer_Smearer.lcsim | 203 +++ .../recon/tracking/kalman/KFOutputDriver.java | 1311 +++++++++++++++++ .../tracking/kalman/KalmanPatRecDriver.java | 13 +- .../mctweaking/RawHitTimeSmearer.java | 231 +++ .../tracking/mctweaking/StripHitSmearer.java | 584 +++----- .../exampleSmearing-0pt008mm.txt | 36 + .../exampleSmearing-0pt01mm.txt | 36 + .../exampleSmearing-0pt0mm.txt | 36 + .../exampleSmearing-0pt1mm.txt | 36 + .../timeSmearing-2ns-4nsL0.txt | 36 + .../timeSmearing-2ns.txt | 36 + .../timeSmearing-5ns.txt | 36 + .../timeSmearing-almost-zero.txt | 36 + 13 files changed, 2205 insertions(+), 425 deletions(-) create mode 100644 steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim create mode 100644 tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java create mode 100644 tracking/src/main/java/org/hps/recon/tracking/mctweaking/RawHitTimeSmearer.java create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt008mm.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt01mm.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt0mm.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt1mm.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns-4nsL0.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-5ns.txt create mode 100644 tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-almost-zero.txt diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim new file mode 100644 index 000000000..c1d734ef0 --- /dev/null +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim @@ -0,0 +1,203 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 12 + false + 10.0 + 15.0 + false + + + + + EcalCalHits EcalClusters EcalClustersCorr FinalStateParticles + UnconstrainedMollerCandidates UnconstrainedMollerVertices UnconstrainedV0Candidates + UnconstrainedV0Vertices TargetConstrainedMollerCandidates TargetConstrainedMollerVertices + TargetConstrainedV0Candidates TargetConstrainedV0Vertices BeamspotConstrainedMollerCandidates + BeamspotConstrainedMollerVertices BeamspotConstrainedV0Candidates BeamspotConstrainedV0Vertices + GBLKinkData GBLKinkDataRelations MatchedToGBLTrackRelations HelicalTrackHits + HelicalTrackHitRelations MatchedTracks GBLTracks MatchedToGBLTrackRelations + RotatedHelicalTrackHits RotatedHelicalTrackHitRelations SVTFittedRawTrackerHits + SVTShapeFitParameters StripClusterer_SiTrackerHitStrip1D TrackData TrackDataRelations + TrackResiduals TrackResidualsRelations OtherElectrons UnconstrainedVcCandidates UnconstrainedVcVertices + KFTrackData KFTrackDataRelations KFGBLStripClusterData KFGBLStripClusterDataRelations SVTShapeFitParameters UnconstrainedVcCandidates_KF UnconstrainedVcVertices_KF + KalmanFullTracks UnconstrainedV0Candidates_KF TargetConstrainedV0Candidates_KF BeamspotConstrainedV0Candidates_KF UnconstrainedV0Vertices_KF TargetConstrainedV0Vertices_KF BeamspotConstrainedV0Vertices_KF + + + + + + CONFIG + + + + + + EcalCalHits + true + 2.4 + + + + + + + WARNING + EcalClusters + + + EcalClusters + EcalClustersCorr + + + + + pass4b-KF-hps2016_L1b_axial_hole.txt pass4b-KF-hps2016_L1b_stereo_hole.txt pass4b-KF-hps2016_L1t_axial_hole.txt pass4b-KF-hps2016_L1t_stereo_hole.txt pass4b-KF-hps2016_L2b_axial_hole.txt pass4b-KF-hps2016_L2b_stereo_hole.txt pass4b-KF-hps2016_L2t_axial_hole.txt pass4b-KF-hps2016_L2t_stereo_hole.txt + false + + + + exampleSmearing-0pt008mm.txt + false + + + + timeSmearing-2ns-4nsL0.txt + false + + + + + + SVTRawTrackerHits + + + Pileup + + true + + 171. + + true + + false + + false + + true + + false + + true + false + + + 24.0 + + + + true + 40 + -4.3 + true + false + + + + EcalClustersCorr + KalmanFullTracks + KalmanFullTracks + FinalStateParticles_KF + UnconstrainedV0Candidates_KF + UnconstrainedV0Vertices_KF + BeamspotConstrainedV0Candidates_KF + BeamspotConstrainedV0Vertices_KF + TargetConstrainedV0Candidates_KF + TargetConstrainedV0Vertices_KF + true + false + -0.224 + 0.275 + -0.08 + 0.060 + -4.3 + 2.8 + 2.8 + 0.0 + 5.0 + 10 + 36.8 + true + true + TrackClusterMatcherMinDistance + + + + + + + ${outputFile}.slcio + + + + + + ${outputFile}.root + false + 0.0 + + KalmanFullTracks + 0.1 + 4.8 + 9999 + true + true + + + + + + diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java new file mode 100644 index 000000000..083c7186c --- /dev/null +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java @@ -0,0 +1,1311 @@ +package org.hps.recon.tracking.kalman; + +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; +import org.hps.util.Pair; + +import java.io.IOException; +import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; +import java.util.logging.Level; +import java.util.logging.Logger; + +import org.hps.recon.tracking.FittedRawTrackerHit; +import org.hps.recon.tracking.ShapeFitParameters; +import org.hps.recon.tracking.CoordinateTransformations; +import org.hps.recon.tracking.TrackStateUtils; +import org.hps.recon.tracking.TrackUtils; +import org.lcsim.detector.ITransform3D; +import org.lcsim.detector.tracker.silicon.ChargeCarrier; +import org.lcsim.detector.tracker.silicon.HpsSiSensor; +import org.lcsim.detector.tracker.silicon.SiSensor; +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.base.BaseTrack; +import org.lcsim.event.TrackState; +import org.lcsim.event.base.BaseTrackState; +import org.lcsim.event.TrackerHit; +import org.lcsim.event.base.BaseRelationalTable; +import org.lcsim.fit.helicaltrack.HelicalTrackCross; +import org.lcsim.geometry.Detector; +import org.lcsim.util.Driver; +import org.lcsim.util.aida.AIDA; +import hep.aida.IManagedObject; +import hep.aida.IBaseHistogram; + +import org.lcsim.fit.helicaltrack.HelixUtils; +import org.lcsim.geometry.FieldMap; + + +// E/p plots +import org.lcsim.event.Cluster; +import org.lcsim.event.ReconstructedParticle; +import org.lcsim.event.TrackState; +//Fiducial cuts on the calorimeter cluster +import org.hps.record.triggerbank.TriggerModule; + +/** + * Make post-KF plots + */ +public class KFOutputDriver extends Driver { + + private AIDA aidaKF; // era public + private String outputPlots = "KalmanTrackingPlots.root"; + private String trackCollectionName = "KalmanTracks"; + private String inputCollectionName = "FinalStateParticles_KF"; + private String trackResidualsRelColName = "KFUnbiasResRelations"; + private Map fittedRawTrackerHitMap = new HashMap(); + private String fittedHitsCollectionName = "SVTFittedRawTrackerHits"; + // private String dataRelationCollection = KFKinkData.DATA_RELATION_COLLECTION; + List _fittedHits; + private List sensors = new ArrayList(); + private double bfield; + public boolean debug = false; + private double chi2Cut = 99999; + private double timeOffset=-40.0; //-40ns for MC ,-55 for data (2016 numbers) + //String kinkFolder = "/kf_kinks/"; + String epullFolder = "/err_pulls/"; + String trkpFolder = "/trk_params/"; + String trkpDetailFolder="/trk_detail/"; + String resFolder="/res/"; + String hitFolder="/hit/"; + String eopFolder = "/EoP/"; + // private boolean b_doKFkinks = false; + private boolean b_doKFresiduals = true; + private boolean b_doDetailPlots = false; + private boolean b_doRawHitPlots = true; + //The field map for extrapolation + private FieldMap bFieldMap; + + //The location of the extrapolation + private double bsZ = 0.; + + //Spacing between top and bottom in the 2D histos + private int mod = 5; + + private double minMom = 1.; + private double maxMom = 6.; + + private double minPhi = -999.9; + private double maxPhi = 999.9; + + private double minTanL = 0.015; + private double maxTanL = 999.9; + + private int nHits = 10; + + private boolean useParticles = true; + + private Pair _trkTimeSigma; + + public void setDebug(boolean val) { + debug = val; + } + + public void setTimeOffset(double val){ + timeOffset=val; + } + + public void setUseParticles(boolean val) { + useParticles = val; + } + /* + public void setDataRelationCollection (String val) { + dataRelationCollection = val; + } + */ + public void setNHits (int val ) { + nHits = val; + } + + public void setMinMom (double val) { + minMom = val; + } + + public void setMaxMom (double val) { + maxMom = val; + } + + public void setMinPhi (double val) { + minPhi = val; + } + + public void setMaxPhi (double val) { + maxPhi = val; + } + + public void setMinTanL (double val) { + minTanL = val; + } + + public void setMaxTanL (double val) { + maxTanL = val; + } + + + //Override the Z of the target. + public void setBsZ (double input) { + bsZ = input; + } + + public void setDoKFresiduals (boolean input) { + b_doKFresiduals = input; + } + + // public void setDoKFkinks (boolean input) { + // b_doKFkinks = input; + // } + + public void setTrackResidualsRelColName (String val) { + trackResidualsRelColName = val; + } + + public void setChi2Cut(double input) { + chi2Cut = input; + } + + public void setOutputPlotsFilename(String fname) { + outputPlots = fname; + } + + public void setTrackCollectionName(String val) { + trackCollectionName=val; + } + + public void setInputCollectionName(String val) { + inputCollectionName=val; + } + + + @Override + protected void detectorChanged(Detector detector) { + if (aidaKF == null) + aidaKF = AIDA.defaultInstance(); + + aidaKF.tree().cd("/"); + + for (HpsSiSensor s : detector.getDetectorElement().findDescendants(HpsSiSensor.class)) { + if (s.getName().startsWith("module_") && s.getName().endsWith("sensor0")) { + sensors.add(s); + } + } + + + Hep3Vector fieldInTracker = TrackUtils.getBField(detector); + this.bfield = Math.abs(fieldInTracker.y()); + + bFieldMap = detector.getFieldMap(); + + if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) { + + // kinkFolder = "/kf_kinks/"; + epullFolder = "/kf_err_pulls/"; + trkpFolder = "/kf_trk_params/"; + trkpDetailFolder = "/kf_trk_detail/"; + resFolder = "/kf_res/"; + hitFolder = "/kf_hit/"; + } + + + + setupPlots(); + setupEoPPlots(); + } + + @Override + public void process(EventHeader event) { + + + + // Track Collection + List tracks = new ArrayList(); + + // Particle Collection + List particles = null; + + // Create a mapping of matched Tracks to corresponding Clusters. + HashMap TrackClusterPairs = new HashMap(); + if(b_doRawHitPlots){ + // Get the list of fitted hits from the event + _fittedHits = event.get(LCRelation.class, fittedHitsCollectionName); + } + int TrackType = 0; + if (!useParticles) { + if (debug) + System.out.println("PF:: DEBUG :: NOT Using particles" + trackCollectionName); + if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) { + TrackType = 1; + } + } + else { + if (debug) + System.out.println("PF:: DEBUG :: Using particles" + inputCollectionName); + if (inputCollectionName.contains("Kalman") || inputCollectionName.contains("KF")) { + + TrackType = 1 ; + } + + } + if (debug) + System.out.println("PF:: DEBUG :: Track Type=" + TrackType); + + + if (!useParticles) + tracks = event.get(Track.class,trackCollectionName); + else { + particles = event.get(ReconstructedParticle.class, inputCollectionName); + for (ReconstructedParticle particle : particles) { + //this requires track cluster match + if (particle.getTracks().isEmpty() || particle.getClusters().isEmpty()) + continue; + Track track = particle.getTracks().get(0); + Cluster cluster = particle.getClusters().get(0); + tracks.add(track); + TrackClusterPairs.put(track,cluster); + } + } + + int nTracks=tracks.size(); + if(debug) + System.out.println(this.getClass()+":: found "+nTracks + " tracks"); + aidaKF.histogram1D(trkpFolder+"nTracks").fill(nTracks); + RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event); + RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event); + + for (Track trk : tracks) { + + if (trk.getChi2() > chi2Cut) + continue; + + if (trk.getTrackerHits().size() < nHits) + continue; + + + if(debug) + System.out.println("Track passed hits d0 = "+trk.getTrackStates().get(0).getD0()); + + Hep3Vector momentum = new BasicHep3Vector(trk.getTrackStates().get(0).getMomentum()); + if (momentum.magnitude() < minMom) + continue; + + if (momentum.magnitude() > maxMom) + continue; + + if(debug) + System.out.println("Track passed momentum"); + + TrackState trackState = trk.getTrackStates().get(0); + if (Math.abs(trackState.getTanLambda()) < minTanL) + continue; + + if (Math.abs(trackState.getTanLambda()) > maxTanL) + continue; + + if (Math.abs(trackState.getPhi()) < minPhi) + continue; + + if (Math.abs(trackState.getPhi()) > maxPhi) + continue; + + if(debug) + System.out.println("Track passed tanLambda"); + + + Map sensorHits = new HashMap(); + + for (TrackerHit hit : trk.getTrackerHits()) { + HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()); + if (sensor != null) { + if(debug){ + System.out.println(this.getClass().getName()+":: inserting hit on sensor = "+sensor.getName()); + } + sensorHits.put(sensor, hit); + } + + if (debug && sensor == null) + System.out.printf("TrackerHit null sensor %s \n", hit.toString()); + } + _trkTimeSigma=getTrackTime(sensorHits); + doBasicKFtrack(trk,sensorHits); + if (b_doKFresiduals) + doKFresiduals(trk, sensorHits,event); + + // if (b_doGBLkinks) + // doGBLkinks(trk,gblKink, sensorNums); + + if (useParticles) + doEoPPlots(trk,TrackClusterPairs.get(trk)); + + + } + } + + private void doEoPPlots(Track track, Cluster cluster) { + + double energy = cluster.getEnergy(); + double[] trk_prms = track.getTrackParameters(); + double tanL = trk_prms[BaseTrack.TANLAMBDA]; + double phi = trk_prms[BaseTrack.PHI]; + TrackState trackState = track.getTrackStates().get(0); + double trackp = new BasicHep3Vector(trackState.getMomentum()).magnitude(); + double eop = energy / trackp; + + String vol = tanL > 0 ? "top" : "bottom"; + + //Charge sign is flipped + String charge = track.getCharge() > 0 ? "ele" : "pos"; + + + aidaKF.histogram1D(eopFolder+"Ecluster_"+vol).fill(energy); + aidaKF.histogram1D(eopFolder+"EoP_"+vol).fill(eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+vol).fill(phi,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+vol).fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+vol).fill(tanL,eop); + Double trackTime = _trkTimeSigma.getFirstElement(); + double trkCluTime=trackTime-cluster.getCalorimeterHits().get(0).getTime()-timeOffset; + aidaKF.histogram1D(trkpFolder+"trk-cluTime_"+charge+"_"+vol).fill(trkCluTime); + aidaKF.histogram1D(trkpFolder+"trk-cluTime_"+vol).fill(trkCluTime); + + + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol).fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol).fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol).fill(phi,eop); + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi").fill(phi,eop); + aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, + phi, + eop); + + + if (TriggerModule.inFiducialRegion(cluster)) { + + aidaKF.histogram1D(eopFolder+"Ecluster_"+vol+"_fid").fill(energy); + aidaKF.histogram1D(eopFolder+"EoP_"+vol+"_fid").fill(eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+vol+"_fid").fill(phi,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+vol+"_fid").fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+vol+"_fid").fill(tanL,eop); + + + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP_"+charge+"_"+vol+"_fid").fill(trackp,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_"+charge+"_"+vol+"_fid").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_"+charge+"_"+vol+"_fid").fill(phi,eop); + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(tanL,eop); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(phi,eop); + aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(tanL, + phi, + eop); + + + + // Cluster positions + + double clusterX = cluster.getPosition()[0]; + double clusterY = cluster.getPosition()[1]; + TrackState ts_ecal = TrackUtils.getTrackStateAtECal(track); + + if(ts_ecal == null){ + return; + } + + double[] ts_ecalPos = ts_ecal.getReferencePoint(); + double trkX = ts_ecalPos[1]; + double trkY = ts_ecalPos[2]; + + aidaKF.histogram1D(eopFolder+"Xcluster_"+vol+"_fid").fill(clusterX); + aidaKF.histogram1D(eopFolder+"Ycluster_"+vol+"_fid").fill(clusterY); + + aidaKF.histogram1D(eopFolder+"trk_clu_resX_"+vol+"_fid").fill(trkX-clusterX); + aidaKF.histogram1D(eopFolder+"trk_clu_resY_"+vol+"_fid").fill(trkY-clusterY); + + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX_"+vol+"_fid").fill(trkX,trkX-clusterX); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY_"+vol+"_fid").fill(trkY,trkX-clusterX); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX_"+vol+"_fid").fill(trkX,trkY-clusterY); + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY_"+vol+"_fid").fill(trkY,trkY-clusterY); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+vol+"_fid").fill(trackp,trkY-clusterY); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+vol+"_fid").fill(trackp,trkX-clusterX); + + aidaKF.histogram2D(eopFolder+"trkY_vs_tanL_"+vol+"_fid").fill(tanL,trkY); + + + aidaKF.histogram1D(eopFolder+"Xcluster_"+charge+"_"+vol+"_fid").fill(clusterX); + aidaKF.histogram1D(eopFolder+"Ycluster_"+charge+"_"+vol+"_fid").fill(clusterY); + + aidaKF.histogram1D(eopFolder+"trk_clu_resX_"+charge+"_"+vol+"_fid").fill(trkX-clusterX); + aidaKF.histogram1D(eopFolder+"trk_clu_resY_"+charge+"_"+vol+"_fid").fill(trkY-clusterY); + + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkX-clusterX); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkX-clusterX); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX_"+charge+"_"+vol+"_fid").fill(trkX,trkY-clusterY); + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY_"+charge+"_"+vol+"_fid").fill(trkY,trkY-clusterY); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkY-clusterY); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP_"+charge+"_"+vol+"_fid").fill(trackp,trkX-clusterX); + + aidaKF.histogram2D(eopFolder+"trkY_vs_tanL_"+charge+"_"+vol+"_fid").fill(tanL,trkY); + + + // + + // As function of incident angle at ECAL, inclusive and in bin of momentum. + + + + + + } + + + } + + + + /* + private void doKFkinks(Track trk, GenericObject kink, Map sensorNums) { + + if (kink == null) { + System.out.println("WARNING::Kink object is null"); + return; + } + + + String vol = "_top"; + int spacing = 0; + if (trk.getTrackStates().get(0).getTanLambda() < 0) { + vol = "_bottom"; + spacing = sensors.size() / 2 + mod; + } + + for (HpsSiSensor sensor : sensorNums.keySet()) { + int index = sensorNums.get(sensor); + double phi = kink.getDoubleVal(index); + float lambda = kink.getFloatVal(index); + + //(2019) For top 0-20, for bottom 25-45 + aidaKF.histogram2D(kinkFolder+"lambda_kink_mod").fill(sensor.getMillepedeId()+spacing,lambda); + aidaKF.profile1D(kinkFolder+"lambda_kink_mod_p").fill(sensor.getMillepedeId()+spacing,lambda); + aidaKF.histogram2D(kinkFolder+"phi_kink_mod").fill(sensor.getMillepedeId()+spacing,phi); + aidaKF.profile1D(kinkFolder+"phi_kink_mod_p").fill(sensor.getMillepedeId()+spacing,phi); + aidaKF.histogram1D(kinkFolder+"lambda_kink_" + sensor.getName()).fill(lambda); + aidaKF.histogram1D(kinkFolder+"phi_kink_" + sensor.getName()).fill(phi); + } + + } + + private void doMTresiduals(Track trk, Map sensorHits) { + TrackState trackState = trk.getTrackStates().get(0); + for (HpsSiSensor sensor : sensorHits.keySet()) { + Hep3Vector extrapPos = TrackStateUtils.getLocationAtSensor(trackState, sensor, bfield); + Hep3Vector hitPos = new BasicHep3Vector(sensorHits.get(sensor).getPosition()); + if (hitPos == null || extrapPos == null) + return; + Hep3Vector diff = VecOp.sub(extrapPos, hitPos); + if (debug) + System.out.printf("MextrapPos %s MhitPos %s \n Mdiff %s ", extrapPos.toString(), hitPos.toString(), diff.toString()); + + ITransform3D trans = sensor.getGeometry().getGlobalToLocal(); + trans.rotate(diff); + + aidaKF.histogram1D(resFolder+"residual_before_KF_" + sensor.getName()).fill(diff.x()); + if (debug) + System.out.printf("MdiffSensor %s \n", diff.toString()); + + } + } + */ + private void FillKFTrackPlot(String str, String isTop, String charge, double val) { + aidaKF.histogram1D(str+isTop).fill(val); + aidaKF.histogram1D(str+isTop+charge).fill(val); + } + + private void FillKFTrackPlot(String str, String isTop, String charge, double valX, double valY) { + aidaKF.histogram2D(str+isTop).fill(valX,valY); + aidaKF.histogram2D(str+isTop+charge).fill(valX,valY); + } + + private void FillKFTrackPlot(String str, String isTop, String charge, double valX, double valY, double valZ) { + aidaKF.histogram3D(str+isTop).fill(valX,valY,valZ); + aidaKF.histogram3D(str+isTop+charge).fill(valX,valY,valZ); + } + + + + private void doBasicKFtrack(Track trk, Map sensorHits) { + + TrackState trackState = trk.getTrackStates().get(0); + + String isTop = "_bottom"; + //if (trk.getTrackerHits().get(0).getPosition()[2] > 0) { + // isTop = "_top"; + //} + + //if (trk.getType()==1 && trk.getTrackerHits().size() < 10) { + // return; + //} + + //List missingHits; + //missingHits = findMissingLayer(trk); + + if (trackState.getTanLambda() > 0) { + isTop = "_top"; + } + + //There is a sign flip in the charge + String charge = "_pos"; + if (trk.getCharge()>0) + charge = "_neg"; + + + //Hep3Vector mom = new BasicHep3Vector(trackState.getMomentum()); + //System.out.println("Track momentum " + mom.toString()); + double trackp = new BasicHep3Vector(trackState.getMomentum()).magnitude(); + + + FillKFTrackPlot(trkpFolder+"d0",isTop,charge,trackState.getD0()); + FillKFTrackPlot(trkpFolder+"z0",isTop,charge,trackState.getZ0()); + FillKFTrackPlot(trkpFolder+"phi",isTop,charge,trackState.getPhi()); + FillKFTrackPlot(trkpFolder+"tanLambda",isTop,charge,trackState.getTanLambda()); + FillKFTrackPlot(trkpFolder+"p",isTop,charge,trackp); + if (trk.getTrackerHits().size()==7) + FillKFTrackPlot(trkpFolder+"p7h",isTop,charge,trackp); + if (trk.getTrackerHits().size()==6) + FillKFTrackPlot(trkpFolder+"p6h",isTop,charge,trackp); + if (trk.getTrackerHits().size()==5) + FillKFTrackPlot(trkpFolder+"p5h",isTop,charge,trackp); + + if (TrackUtils.isHoleTrack(trk)) + FillKFTrackPlot(trkpFolder+"p_hole",isTop,charge,trackp); + else + FillKFTrackPlot(trkpFolder+"p_slot",isTop,charge,trackp); + Double trackTime = _trkTimeSigma.getFirstElement(); + Double trackTimeSD = _trkTimeSigma.getSecondElement(); + + //fill track time and standard dev + FillKFTrackPlot(trkpFolder+"trkTime",isTop,charge,trackTime); + FillKFTrackPlot(trkpFolder+"trkTimeSD",isTop,charge,trackTimeSD); + + //Momentum maps + FillKFTrackPlot(trkpFolder+"p_vs_phi",isTop,charge,trackState.getPhi(),trackp); + FillKFTrackPlot(trkpFolder+"p_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackp); + FillKFTrackPlot(trkpFolder+"p_vs_phi_tanLambda",isTop,charge,trackState.getPhi(),trackState.getTanLambda(),trackp); + + double tanLambda = trackState.getTanLambda(); + double cosLambda = 1. / (Math.sqrt(1+tanLambda*tanLambda)); + + FillKFTrackPlot(trkpFolder+"pT_vs_phi",isTop,charge,trackState.getPhi(),trackp*cosLambda); + FillKFTrackPlot(trkpFolder+"pT_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackp*cosLambda); + + + //if (trk.getTrackerHits().size()==6) + // FillKFTrackPlot(trkpFolder+"p_Missing1Hit",isTop,charge,missingHits.get(0),trackp); + + //if (missingHits.size()==1 && missingHits.get(0)==7) + // FillKFTrackPlot(trkpFolder+"p_MissingLastLayer",isTop,charge,trackp); + + FillKFTrackPlot(trkpFolder+"Chi2",isTop,charge,trk.getChi2()); + FillKFTrackPlot(trkpFolder+"Chi2oNDF",isTop,charge,trk.getChi2() / trk.getNDF()); + FillKFTrackPlot(trkpFolder+"Chi2_vs_p",isTop,charge,trackp,trk.getChi2()); + + int nhits = trk.getTrackerHits().size(); + + aidaKF.histogram1D(trkpFolder+"nHits" + isTop).fill(nhits); + aidaKF.histogram1D(trkpFolder+"nHits" + isTop+charge).fill(nhits); + + Hep3Vector beamspot = CoordinateTransformations.transformVectorToDetector(TrackUtils.extrapolateHelixToXPlane(trackState, 0)); + if (debug) + System.out.printf("beamspot %s transformed \n", beamspot.toString()); + FillKFTrackPlot(trkpFolder+"trk_extr_or_x",isTop,charge,beamspot.x()); + FillKFTrackPlot(trkpFolder+"trk_extr_or_y",isTop,charge,beamspot.y()); + + //Extrapolation to assumed tgt pos - helix + Hep3Vector trkTgt = CoordinateTransformations.transformVectorToDetector(TrackUtils.extrapolateHelixToXPlane(trackState,bsZ)); + FillKFTrackPlot(trkpFolder+"trk_extr_bs_x",isTop,charge,trkTgt.x()); + FillKFTrackPlot(trkpFolder+"trk_extr_bs_y",isTop,charge,trkTgt.y()); + + //Transform z to the beamspot plane + //Get the PathToPlane + + BaseTrackState ts_bs = TrackUtils.getTrackExtrapAtVtxSurfRK(trackState,bFieldMap,0.,bsZ); + + + //Get the track parameters wrt the beamline using helix + double [] beamLine = new double [] {bsZ,0}; + double [] helixParametersAtBS = TrackUtils.getParametersAtNewRefPoint(beamLine, trackState); + + + FillKFTrackPlot(trkpFolder+"trk_extr_bs_x_rk",isTop,charge,ts_bs.getReferencePoint()[1]); + FillKFTrackPlot(trkpFolder+"trk_extr_bs_y_rk",isTop,charge,ts_bs.getReferencePoint()[2]); + + //Ill defined - should be defined wrt bsX and bsY + FillKFTrackPlot(trkpFolder+"d0_vs_bs_rk",isTop,charge,ts_bs.getD0()); + FillKFTrackPlot(trkpFolder+"d0_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.D0]); + + double s = HelixUtils.PathToXPlane(TrackUtils.getHTF(trackState),bsZ,0.,0).get(0); + FillKFTrackPlot(trkpFolder+"z0_vs_bs",isTop,charge,trackState.getZ0() + s*trackState.getTanLambda()); + FillKFTrackPlot(trkpFolder+"z0_vs_bs_rk",isTop,charge,ts_bs.getZ0()); + FillKFTrackPlot(trkpFolder+"z0_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.Z0]); + + + FillKFTrackPlot(trkpFolder+"phi_vs_bs_extrap",isTop,charge,helixParametersAtBS[BaseTrack.PHI]); + + //TH2D - Filling + FillKFTrackPlot(trkpFolder+"d0_vs_phi",isTop,charge,trackState.getPhi(),trackState.getD0()); + FillKFTrackPlot(trkpFolder+"d0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getD0()); + FillKFTrackPlot(trkpFolder+"d0_vs_p",isTop,charge,trackp,trackState.getD0()); + + //Ill defined - should be defined wrt bsX and bsY + FillKFTrackPlot(trkpFolder+"d0bs_vs_p",isTop,charge,trackp,helixParametersAtBS[BaseTrack.D0]); + + FillKFTrackPlot(trkpFolder+"z0_vs_p",isTop,charge,trackp,trackState.getZ0()); + FillKFTrackPlot(trkpFolder+"z0bs_vs_p",isTop,charge,trackp,ts_bs.getZ0()); + + //Interesting plot to get a sense where z-vtx is. + //If z0 is referenced to the right BS z location, the slope of vs tanLambda is 0 + FillKFTrackPlot(trkpFolder+"z0_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackState.getZ0()); + FillKFTrackPlot(trkpFolder+"z0bs_vs_tanLambda",isTop,charge,trackState.getTanLambda(),ts_bs.getZ0()); + + + if(b_doRawHitPlots){ + + // Map the fitted hits to their corresponding raw hits + this.mapFittedRawHits(_fittedHits); + + for(TrackerHit tkh: trk.getTrackerHits()){ + List rawhits = tkh.getRawHits(); + for(RawTrackerHit rth: rawhits){ + //need the rth->fited + HpsSiSensor sensor = (HpsSiSensor) rth.getDetectorElement(); + double t0 = FittedRawTrackerHit.getT0(getFittedHit(rth)); + double amplitude = FittedRawTrackerHit.getAmp(getFittedHit(rth)); + double chi2Prob = ShapeFitParameters.getChiProb(FittedRawTrackerHit.getShapeFitParameters(getFittedHit(rth))); + aidaKF.histogram1D(hitFolder+"raw_hit_t0_"+sensor.getName()).fill(t0); + aidaKF.histogram1D(hitFolder+"raw_hit_amplitude_"+sensor.getName()).fill(amplitude); + aidaKF.histogram1D(hitFolder+"raw_hit_chisq_"+sensor.getName()).fill(chi2Prob); + + } + } + } + if (b_doDetailPlots) { + int ibins = 15; + double start= -12; + double end = -5; + double step = (end-start) / (double)ibins; + + for (int ibin = 0; ibin sensorHits, EventHeader event) { + + Map sensorMPIDs = new HashMap(); + + for (HpsSiSensor sensor : sensorHits.keySet()) { + //Also fill here the sensorMPIDs map + if(debug){ + System.out.println(this.getClass().getName()+":: mapping "+sensor.getMillepedeId()+" to " + sensor.getName()); + } + sensorMPIDs.put(sensor.getMillepedeId(),sensor); + ITransform3D trans = sensor.getGeometry().getGlobalToLocal(); + + // position of hit (track crossing the sensor before kf extrapolation) + // the hit information available on each sensor is meaningful only along the measurement direction, + // Hep3Vector hitPos = new BasicHep3Vector(sensorHits.get(sensor).getPosition()); + // instead: extract the information of the hit of the track at the sensor position before kf + TrackState trackState = trk.getTrackStates().get(0); + Hep3Vector hitTrackPos = TrackStateUtils.getLocationAtSensor(trackState, sensor, bfield); + + if (hitTrackPos == null) { + if (debug) { + System.out.printf(this.getClass().getName()+"::doKFresiduals:: hitTrackPos is null to sensor %s\n", sensor.toString()); + } + continue; + } + + Hep3Vector hitTrackPosSensor = new BasicHep3Vector(hitTrackPos.v()); + trans.transform(hitTrackPosSensor); + // after the transformation x and y in the sensor frame are reversed + // This plot is ill defined. + + aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitTrackPosSensor.y(), hitTrackPosSensor.x()); + //aidaKF.histogram2D("hit_u_vs_v_sensor_frame_" + sensor.getName()).fill(hitPos.y(), hitPos.x()); + //aidaKF.histogram2D("hit y vs x lab-frame " + sensor.getName()).fill(hitPos.y(), hitPos.x()); + + + // position predicted on track after KF + Hep3Vector extrapPos = null; + Hep3Vector extrapPosSensor = null; + extrapPos = TrackUtils.extrapolateTrackPositionToSensor(trk, sensor, sensors, bfield); + if (extrapPos == null) + return; + extrapPosSensor = new BasicHep3Vector(extrapPos.v()); + trans.transform(extrapPosSensor); + //aidaKF.histogram2D("residual after KF vs u predicted " + sensor.getName()).fill(extrapPosSensor.x(), res); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + // select track charge + if(trk.getCharge()>0) { + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + }else if(trk.getCharge()<0) { + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName()).fill(extrapPosSensor.y(), extrapPosSensor.x()); + } + + // post-KF residual + Hep3Vector hitPos = new BasicHep3Vector(sensorHits.get(sensor).getPosition()); + Hep3Vector hitPosSensor = new BasicHep3Vector(hitPos.v()); + trans.transform(hitPosSensor); + Hep3Vector resSensor = VecOp.sub(hitPosSensor, extrapPosSensor); + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName()).fill(extrapPosSensor.y(), resSensor.x()); + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName()).fill(hitPosSensor.x(), resSensor.x()); + aidaKF.histogram1D(resFolder+"residual_after_KF_" + sensor.getName()).fill(resSensor.x()); + + + + if (debug) { + System.out.printf("hitPos %s hitPosSensor %s \n", hitPos.toString(), hitPosSensor.toString()); + System.out.printf("resSensor %s \n", resSensor.toString()); + System.out.printf("extrapPos %s extrapPosSensor %s \n", extrapPos.toString(), extrapPosSensor.toString()); + ITransform3D electrodes_to_global = sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getLocalToGlobal(); + Hep3Vector measuredCoordinate = sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getMeasuredCoordinate(0); + Hep3Vector unmeasuredCoordinate = sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getUnmeasuredCoordinate(0); + System.out.printf("unMeasCoordOrig %s MeasCoordOrig %s \n", unmeasuredCoordinate.toString(), measuredCoordinate.toString()); + measuredCoordinate = VecOp.mult(VecOp.mult(CoordinateTransformations.getMatrix(), electrodes_to_global.getRotation().getRotationMatrix()), measuredCoordinate); + unmeasuredCoordinate = VecOp.mult(VecOp.mult(CoordinateTransformations.getMatrix(), electrodes_to_global.getRotation().getRotationMatrix()), unmeasuredCoordinate); + Hep3Vector testX = trans.inverse().rotated(new BasicHep3Vector(1, 0, 0)); + Hep3Vector testY = trans.inverse().rotated(new BasicHep3Vector(0, 1, 0)); + Hep3Vector testZ = trans.inverse().rotated(new BasicHep3Vector(0, 0, 1)); + System.out.printf("unMeasCoord %s MeasCoord %s \n transX %s transY %s transZ %s \n", unmeasuredCoordinate.toString(), measuredCoordinate.toString(), testX.toString(), testY.toString(), testZ.toString()); + } + }//loop on sensor hits + + Double trackTime = _trkTimeSigma.getFirstElement(); + Double trackTimeSD = _trkTimeSigma.getSecondElement(); + /* + trackTime /= (float)sensorHits.size(); + + for (HpsSiSensor sensor : sensorHits.keySet()) { + trackTimeSD += Math.pow(trackTime - sensorHits.get(sensor).getTime(),2); + } + + trackTimeSD = Math.sqrt(trackTimeSD / ((float) sensorHits.size() - 1.)); + */ + + + RelationalTable trackResidualsTable = null; + if (event.hasCollection(LCRelation.class, trackResidualsRelColName)) { + trackResidualsTable = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED); + List trackresRelation = event.get(LCRelation.class, trackResidualsRelColName); + for (LCRelation relation : trackresRelation) { + if (relation != null && relation.getFrom() != null && relation.getTo() != null) { + trackResidualsTable.add(relation.getFrom(), relation.getTo()); + } + } + if (debug) + System.out.println("Loaded track Residuals Table"); + } else { + if (debug) { + System.out.println("null TrackResidualsKF Data Relations."); + } + //Failed finding TrackResidualsKF + return; + } + + GenericObject trackRes = (GenericObject) trackResidualsTable.from(trk); + if (trackRes == null) { + if (debug) + System.out.println("null TrackResidualsKF Data."); + return; + } + + int nres = (trackRes.getNInt()-1); + if(debug){ + System.out.println(this.getClass().getName()+":: number entries in trackRes = "+nres); + } + String vol = "_top"; + if (trk.getTrackStates().get(0).getTanLambda() < 0) + vol = "_bottom"; + // get the unbias + for (int i_hit =0; i_hit < nres ; i_hit+=1) { + if (trackRes.getIntVal(i_hit)!=-999) { + if(debug){ + System.out.println(this.getClass().getName()+":: getting residual for ihit = "+i_hit+" trackResValue = "+ trackRes.getIntVal(i_hit)); + } + //Measured hit + HpsSiSensor hps_sensor = sensorMPIDs.get(trackRes.getIntVal(i_hit)); + if(debug){ + System.out.println(this.getClass().getName()+":: sensor for this hit = "+hps_sensor.getName()); + } + Hep3Vector hitPosG = new BasicHep3Vector(sensorHits.get(hps_sensor).getPosition()); + Hep3Vector hitPosSensorG = new BasicHep3Vector(hitPosG.v()); + ITransform3D g2l = hps_sensor.getGeometry().getGlobalToLocal(); + g2l.transform(hitPosSensorG); + String sensorName = (sensorMPIDs.get(trackRes.getIntVal(i_hit))).getName(); + + //Predicted hit + Hep3Vector extrapPos = null; + Hep3Vector extrapPosSensor = null; + extrapPos = TrackUtils.extrapolateTrackPositionToSensor(trk, hps_sensor, sensors, bfield); + if (extrapPos == null) + continue; + extrapPosSensor = new BasicHep3Vector(extrapPos.v()); + g2l.transform(extrapPosSensor); + + if (debug) { + System.out.printf("NHits %d MPID sensor:%d %s %d\n", nres,trackRes.getIntVal(i_hit), sensorName,i_hit); + System.out.printf("Track uresiduals: %s %.5f %.5f\n",sensorName, trackRes.getDoubleVal(i_hit),trackRes.getFloatVal(i_hit)); + } + + //General residuals Per volume + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol).fill(trackRes.getDoubleVal(i_hit)); + + if (trackRes.getIntVal(i_hit) < 9) + //L1L4 + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol+"_L1L4").fill(trackRes.getDoubleVal(i_hit)); + else + //L5L7 + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol+"_L5L7").fill(trackRes.getDoubleVal(i_hit)); + + + //Top go from 0 to 20, bottom go from 25 to 45 + int spacing = 0; + if (vol == "_bottom") + spacing = sensors.size()/2 + mod; + + aidaKF.histogram2D(resFolder+"uresidual_KF_mod").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); + aidaKF.profile1D(resFolder+"uresidual_KF_mod_p").fill(trackRes.getIntVal(i_hit)+spacing,trackRes.getDoubleVal(i_hit)); + aidaKF.histogram1D(resFolder+"uresidual_KF_" + sensorName).fill(trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensorName).fill(hitPosSensorG.x(),trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensorName).fill(extrapPosSensor.y(),trackRes.getDoubleVal(i_hit)); + aidaKF.histogram1D(epullFolder+"ureserror_KF_" + sensorName).fill(trackRes.getFloatVal(i_hit)); + aidaKF.histogram1D(epullFolder+"ures_pull_KF_" + sensorName).fill(trackRes.getDoubleVal(i_hit) / Math.sqrt(trackRes.getFloatVal(i_hit))); + + //Get the hit time + double hitTime = sensorHits.get(hps_sensor).getTime(); + + //Get the track time (it's the average of hits-on-track time) + + double dT_hit_track = hitTime - trackTime; + double dT_hit_sigma = (hitTime - trackTime) / trackTimeSD; + + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_"+sensorName).fill(dT_hit_track,trackRes.getDoubleVal(i_hit)); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_"+sensorName).fill(dT_hit_sigma,trackRes.getDoubleVal(i_hit)); + + + + + } + else { + if (debug){ + System.out.printf("Track refit failed? No biased residual for %d\n", i_hit); + } + } + } + }//doKFresiduals + + private List findMissingLayer(Track trk) { + + List layers = new ArrayList(); + layers.add(1); + layers.add(2); + layers.add(3); + layers.add(4); + layers.add(5); + layers.add(6); + layers.add(7); + + List LayersOnTrack = new ArrayList(); + List missingHits = new ArrayList(); + + for (TrackerHit hit : trk.getTrackerHits()) { + int stripLayer = ((HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement()).getLayerNumber(); + // int hpslayer = (stripLayer + 1 ) / 2; + LayersOnTrack.add(stripLayer); + } + for (Integer layer : layers) { + if (!LayersOnTrack.contains(layer)) + missingHits.add(layer); + } + return missingHits; + } + + private void setupEoPPlots() { + + List volumes = new ArrayList(); + volumes.add("_top"); + volumes.add("_bottom"); + + List charges = new ArrayList(); + charges.add(""); + charges.add("_ele"); + charges.add("_pos"); + + for (String vol : volumes) { + + aidaKF.histogram1D(eopFolder+"Ecluster"+vol,200,0,6); + aidaKF.histogram1D(eopFolder+"EoP"+vol,200,0,2); + + double lmin = 0.; + double lmax = 0.08; + if (vol == "_bot") { + lmin = -0.08; + lmax = 0.; + } + + for (String charge : charges) { + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP"+charge+vol,200,0,6,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda"+charge+vol,200,lmin,lmax,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi"+charge+vol,200,-0.2,0.2,200,0,2); + } + + aidaKF.histogram1D(eopFolder+"Ecluster"+vol+"_fid",200,0,5); + aidaKF.histogram1D(eopFolder+"EoP"+vol+"_fid",200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP"+vol+"_fid",200,0,6,200,0,2); + + + double cxrange = 20; + double cyrange = 20; + double ecalX = 400; + + aidaKF.histogram1D(eopFolder+"Xcluster"+vol+"_fid",200,-ecalX,ecalX); + aidaKF.histogram1D(eopFolder+"Ycluster"+vol+"_fid",200,-ecalX,ecalX); + aidaKF.histogram1D(eopFolder+"trk_clu_resX"+vol+"_fid",200,-cxrange,cxrange); + aidaKF.histogram1D(eopFolder+"trk_clu_resY"+vol+"_fid",200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX"+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY"+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX"+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY"+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP"+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP"+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trkY_vs_tanL"+vol+"_fid",200,-0.2,0.2,200,-100,100); + + + for (String charge : charges) { + + //put the trk-cluster time in trkpFolder + aidaKF.histogram1D(trkpFolder+"trk-cluTime"+charge+vol,100,-20,20); + + aidaKF.histogram2D(eopFolder+"EoP_vs_trackP"+charge+vol+"_fid",200,0,6,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda"+charge+vol+"_fid",200,0.01,0.08,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi"+charge+vol+"_fid",200,-0.2,0.2,200,0,2); + + + + aidaKF.histogram1D(eopFolder+"Xcluster"+charge+vol+"_fid",200,-ecalX,ecalX); + aidaKF.histogram1D(eopFolder+"Ycluster"+charge+vol+"_fid",200,-ecalX,ecalX); + aidaKF.histogram1D(eopFolder+"trk_clu_resX"+charge+vol+"_fid",200,-cxrange,cxrange); + aidaKF.histogram1D(eopFolder+"trk_clu_resY"+charge+vol+"_fid",200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsX"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vsY"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cxrange,cxrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsX"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vsY"+charge+vol+"_fid",200,-ecalX,ecalX,200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trk_clu_resY_vstrkP"+charge+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + aidaKF.histogram2D(eopFolder+"trk_clu_resX_vstrkP"+charge+vol+"_fid",100,0.,5,200,-cyrange,cyrange); + + aidaKF.histogram2D(eopFolder+"trkY_vs_tanL"+charge+vol+"_fid",200,-0.2,0.2,200,-100,100); + + + + + } + } + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda",200,-0.1,0.1,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi",200,-0.2,0.2,200,0,2); + aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + + aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid",200,-0.1,0.1,200,0,2); + aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid",200,-0.2,0.2,200,0,2); + aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + + } + + + private void setupPlots() { + + + double xmax = 0.25; + double kxmax = 0.001; + + int nbins = 250; + List volumes = new ArrayList(); + volumes.add("_top"); + volumes.add("_bottom"); + int mod_2dplot_bins = sensors.size()+mod*2; + + for (String vol : volumes) { + // aidaKF.histogram1D(resFolder+"bresidual_KF"+vol,nbins, -xmax, xmax); + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol,nbins, -xmax, xmax); + // aidaKF.histogram1D(resFolder+"bresidual_KF"+vol+"_L1L4",nbins,-xmax,xmax); + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol+"_L1L4",nbins,-xmax,xmax); + // aidaKF.histogram1D(resFolder+"bresidual_KF"+vol+"_L5L7",nbins,-xmax,xmax); + aidaKF.histogram1D(resFolder+"uresidual_KF"+vol+"_L5L7",nbins,-xmax,xmax); + + } + + //res/kinks TH2D + //5 empty bins to distinguish between top and bottom + + // aidaKF.histogram2D(resFolder+"bresidual_KF_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5, nbins, -xmax,xmax); + // aidaKF.profile1D(resFolder+"bresidual_KF_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); + aidaKF.histogram2D(resFolder+"uresidual_KF_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5, 400, -0.4,0.4); + aidaKF.profile1D(resFolder+"uresidual_KF_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); + + + //Hits vs channel + /* + int nch = 400; + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL1b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL2b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL3b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL4b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL5b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL6b",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL7b",nch,0,nch,nch,0,nch); + + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL1t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL2t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL3t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL4t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL5t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL6t",nch,0,nch,nch,0,nch); + aidaKF.histogram2D(resFolder+"Axial_vs_Stereo_channel_moduleL7t",nch,0,nch,nch,0,nch); + */ + + + for (SiSensor sensor : sensors) { + + HpsSiSensor sens = (HpsSiSensor) sensor.getGeometry().getDetectorElement(); + xmax = 0.5; + nbins = 250; + int l = (sens.getLayerNumber() + 1) / 2; + if (l > 1) xmax = 0.05 + (l - 1) * 0.08; + // aidaKF.histogram1D(resFolder+"residual_before_KF_" + sensor.getName(), nbins, -xmax, xmax); + + xmax = 0.250; + + if (l >= 6) + xmax = 0.250; + aidaKF.histogram1D(resFolder+"residual_after_KF_" + sensor.getName(), nbins, -xmax, xmax); + // aidaKF.histogram1D(resFolder+"bresidual_KF_" + sensor.getName(), nbins, -xmax, xmax); + aidaKF.histogram1D(resFolder+"uresidual_KF_" + sensor.getName(), nbins, -xmax, xmax); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_u_hit_" + sensor.getName(),100,-20.0,20.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_v_pred_" + sensor.getName(),300,-60.0,60.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dT_hit_" + sensor.getName(),100,-10.0,10.0,100,-0.1,0.1); + aidaKF.histogram2D(resFolder+"uresidual_KF_vs_dTs_hit_" + sensor.getName(),100,-5.0,5.0,100,-0.1,0.1); + + + // aidaKF.histogram1D(epullFolder+"breserror_KF_" + sensor.getName(), nbins, 0.0, 0.1); + aidaKF.histogram1D(epullFolder+"ureserror_KF_" + sensor.getName(), nbins, 0.0, 0.2); + // aidaKF.histogram1D(epullFolder+"bres_pull_KF_" + sensor.getName(), nbins, -5, 5); + aidaKF.histogram1D(epullFolder+"ures_pull_KF_" + sensor.getName(), nbins, -5, 5); + + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_u_hit_" + sensor.getName(), 100, -20.0, 20.0, 100, -0.04, 0.04); + aidaKF.histogram2D(resFolder+"residual_after_KF_vs_v_predicted_" + sensor.getName(), 100, -55.0, 55.0, 100, -0.04, 0.04); + aidaKF.histogram2D(hitFolder+"hit_u_vs_v_sensor_frame_" + sensor.getName(), 300, -60.0, 60.0, 300, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_pos_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + aidaKF.histogram2D(hitFolder+"predicted_u_vs_v_neg_sensor_frame_" + sensor.getName(), 100, -60, 60, 100, -25, 25); + + aidaKF.histogram1D(hitFolder+"raw_hit_t0_"+sensor.getName(),200, -100, 100.0); + aidaKF.histogram1D(hitFolder+"raw_hit_amplitude_"+sensor.getName(),200, 0.0, 4000.0); + aidaKF.histogram1D(hitFolder+"raw_hit_chisq_"+sensor.getName(),200, 0.0, 2.0); + + + + xmax = 0.0006; + if(l==1){ + xmax = 0.0002; + }else if(l==2){ + xmax = 0.0005; + }else if(l==3 || l==4){ + xmax = 0.0006; + }else if(l >= 5) { + if (sens.isBottomLayer() && sens.isAxial()) + xmax = 0.001; + if (sens.isTopLayer() && !sens.isAxial()) + xmax = 0.001; + } + // aidaKF.histogram1D(kinkFolder+"lambda_kink_" + sensor.getName(), 250, -xmax, xmax); + //aidaKF.histogram1D(kinkFolder+"phi_kink_" + sensor.getName(), 250, -xmax, xmax); + } + /* + aidaKF.histogram2D(kinkFolder+"lambda_kink_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5,nbins,-0.001,0.001); + aidaKF.profile1D(kinkFolder+"lambda_kink_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); + aidaKF.histogram2D(kinkFolder+"phi_kink_mod",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5 ,nbins,-0.001,0.001); + aidaKF.profile1D(kinkFolder+"phi_kink_mod_p",mod_2dplot_bins,-0.5,mod_2dplot_bins-0.5); + */ + List charges = new ArrayList(); + charges.add(""); + charges.add("_pos"); + charges.add("_neg"); + + int nbins_t = 200; + + //For momentum + int nbins_p = 150; + double pmax = 4.; + + double z0max = 1; + double d0max = 5; + double z0bsmax = 0.2; + + aidaKF.histogram1D(trkpFolder+"nTracks",15,0,15); + for (String vol : volumes) { + for (String charge : charges) { + + + //TH1Ds + aidaKF.histogram1D(trkpFolder+"d0"+vol+charge,nbins_t,-5.0,5.0); + aidaKF.histogram1D(trkpFolder+"z0"+vol+charge,nbins_t,-1.3,1.3); + aidaKF.histogram1D(trkpFolder+"phi"+vol+charge,nbins_t,-0.06,0.06); + aidaKF.histogram1D(trkpFolder+"tanLambda"+vol+charge,nbins_t,-0.2,0.2); + aidaKF.histogram1D(trkpFolder+"trkTime"+vol+charge,nbins_t,-20,20); + aidaKF.histogram1D(trkpFolder+"trkTimeSD"+vol+charge,nbins_t,0,10); + + aidaKF.histogram1D(trkpFolder+"p"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p7h"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p6h"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p5h"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p_MissingLastLayer"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p_hole"+vol+charge,nbins_p,0.,pmax); + aidaKF.histogram1D(trkpFolder+"p_slot"+vol+charge,nbins_p,0.,pmax); + + aidaKF.histogram1D(trkpFolder+"Chi2"+vol+charge,nbins_t*2,0,200); + aidaKF.histogram1D(trkpFolder+"Chi2oNDF"+vol+charge,nbins_t*2,0,50); + aidaKF.histogram1D(trkpFolder+"nHits"+vol+charge,15,0,15); + aidaKF.histogram1D(trkpFolder+"trk_extr_or_x"+vol+charge,nbins_t,-3,3); + aidaKF.histogram1D(trkpFolder+"trk_extr_or_y"+vol+charge,nbins_t,-3,3); + aidaKF.histogram1D(trkpFolder+"trk_extr_bs_x"+vol+charge, 2*nbins_t, -5, 5); + aidaKF.histogram1D(trkpFolder+"trk_extr_bs_y"+vol+charge, 2*nbins_t, -5, 5); + aidaKF.histogram1D(trkpFolder+"trk_extr_bs_x_rk"+vol+charge, 2*nbins_t, -5, 5); + aidaKF.histogram1D(trkpFolder+"trk_extr_bs_y_rk"+vol+charge, 2*nbins_t, -3, 3); + aidaKF.histogram1D(trkpFolder+"d0_vs_bs_rk"+vol+charge, 2*nbins_t, -5, 5); + aidaKF.histogram1D(trkpFolder+"d0_vs_bs_extrap"+vol+charge, 2*nbins_t, -5, 5); + aidaKF.histogram1D(trkpFolder+"z0_vs_bs_rk"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); + aidaKF.histogram1D(trkpFolder+"z0_vs_bs_extrap"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); + aidaKF.histogram1D(trkpFolder+"z0_vs_bs"+vol+charge, 2*nbins_t, -z0bsmax, z0bsmax); + aidaKF.histogram1D(trkpFolder+"phi_vs_bs_extrap"+vol+charge,2*nbins_t, -0.06,0.06); + + + //TH2Ds + + aidaKF.histogram2D(trkpFolder+"d0_vs_phi"+vol+charge,nbins_t,-0.3,0.3,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"Chi2_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t*2,0,200); + //aidaKF.histogram2D("d0_vs_phi_bs"+vol+charge,nbins_t,-5.0,5.0,nbins_t,-0.3,0.3); + aidaKF.histogram2D(trkpFolder+"d0_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"d0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"d0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"z0_vs_p"+vol+charge, nbins_p,0.0,pmax,nbins_t,-5.0,5.0); + aidaKF.histogram2D(trkpFolder+"z0bs_vs_p"+vol+charge,nbins_p,0.0,pmax,nbins_t,-z0bsmax,z0bsmax); + aidaKF.histogram2D(trkpFolder+"z0_vs_tanLambda"+vol+charge, nbins_t,-0.1,0.1,nbins_t,-z0max,z0max); + aidaKF.histogram2D(trkpFolder+"z0bs_vs_tanLambda"+vol+charge,nbins_t,-0.1,0.1,nbins_t,-z0bsmax,z0bsmax); + + aidaKF.histogram2D(trkpFolder+"p_Missing1Hit"+vol+charge,8,0,8,nbins_p,0.0,pmax); + aidaKF.histogram2D(trkpFolder+"p_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaKF.histogram2D(trkpFolder+"p_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + aidaKF.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); + + aidaKF.histogram2D(trkpFolder+"pT_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); + aidaKF.histogram2D(trkpFolder+"pT_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); + + + + if (b_doDetailPlots) { + //TH2Ds - detail + int ibins = 15; + double start= -12; + double end = -5; + double step = (end-start) / (double)ibins; + for (int ibin = 0; ibin getTrackTime(Map sensorHits){ + double trackTime = 0.; + double trackTimeSD = 0.; + for (HpsSiSensor sensor : sensorHits.keySet()) { + trackTime += sensorHits.get(sensor).getTime(); + } + trackTime /= (float)sensorHits.size(); + + for (HpsSiSensor sensor : sensorHits.keySet()) { + trackTimeSD += Math.pow(trackTime - sensorHits.get(sensor).getTime(),2); + } + + trackTimeSD = Math.sqrt(trackTimeSD / ((float) sensorHits.size() - 1.)); + return new Pair(trackTime, trackTimeSD); + } + + public void endOfData() { + if (outputPlots != null) { + try { + aidaKF.saveAs(outputPlots); + + /* + // remove all KF histograms from heap after they have been written on output file + String[] type = aidaKF.tree().listObjectNames("/",true); + for (int i=0; i allHits) { + + // Clear the fitted raw hit map of old values + fittedRawTrackerHitMap.clear(); + + // Loop through all fitted hits and map them to their corresponding raw hits + for (LCRelation fittedHit : allHits) { + fittedRawTrackerHitMap.put(FittedRawTrackerHit.getRawTrackerHit(fittedHit), fittedHit); + } + } + + + private LCRelation getFittedHit(RawTrackerHit rawHit) { + return fittedRawTrackerHitMap.get(rawHit); + } +} diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java index e1e85b2d3..6121d5e67 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java @@ -510,21 +510,24 @@ private ArrayList[] prepareTrackCollections(EventHeader event, List res_and_sigma = kTk.unbiasedResidual(ilay); + //...loop over clusters and save millipedID to residuals + for (GBLStripClusterData clstr : clstrs) { + Pair res_and_sigma = kTk.unbiasedResidualMillipede(clstr.getId()); if (res_and_sigma.getSecondElement() > -1.) { - layers.add(ilay); + layers.add(clstr.getId()); residuals.add(res_and_sigma.getFirstElement()); sigmas.add(res_and_sigma.getSecondElement().floatValue()); } - Pair inter_and_sigma = kTk.unbiasedIntersect(ilay, true); + }//Loop on clusters + for (int ilay = 0; ilay<14; ilay++) { + Pair inter_and_sigma = kTk.unbiasedIntersect(ilay, true); layersInt.add(ilay); intersect.add(inter_and_sigma.getFirstElement()[uindex]); intersect.add(inter_and_sigma.getFirstElement()[vindex]); intersect.add(inter_and_sigma.getFirstElement()[windex]); sigmasInt.add(inter_and_sigma.getSecondElement().floatValue()); }//Loop on layers - + //Add the Track Data TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY, (float) svtCenter_bFieldY); trackDataCollection.add(KFtrackData); diff --git a/tracking/src/main/java/org/hps/recon/tracking/mctweaking/RawHitTimeSmearer.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/RawHitTimeSmearer.java new file mode 100644 index 000000000..c022df664 --- /dev/null +++ b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/RawHitTimeSmearer.java @@ -0,0 +1,231 @@ +package org.hps.recon.tracking; + +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStream; +import java.io.InputStreamReader; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.Random; +import java.util.logging.Level; +import java.util.logging.Logger; +import java.util.regex.Matcher; +import java.util.regex.Pattern; +import org.lcsim.detector.tracker.silicon.ChargeCarrier; +import org.lcsim.detector.tracker.silicon.HpsSiSensor; +import org.lcsim.detector.tracker.silicon.SiSensorElectrodes; +import org.lcsim.event.EventHeader; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.TrackerHit; +import org.lcsim.event.base.BaseTrackerHit; +import org.lcsim.geometry.Detector; +import org.lcsim.lcio.LCIOUtil; +import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D; +import org.lcsim.util.Driver; +import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType; +import hep.physics.matrix.SymmetricMatrix; +import org.lcsim.event.LCRelation; + +/** + * + * @author mgraham created 9/20/24 + * This driver will randomize the position of strip clusters + * within a specified gaussian sigma. + * 11/17/24 + * add in cluster time smearing as well. + * + * Strip collection: + * "StripClusterer_SiTrackerHitStrip1D" (default) + * Important Settables: + * smear: the sigma of smearing gaussian + * + * mg...right now, I smear all hits in all layers by the same amount + * independent of nHits. + */ +public class RawHitTimeSmearer extends Driver { + + + String smearTimeFile="foobar"; + //IMPORTANT...the layer, top/bottom/stereo/axial/slot/hole are derived from these names!!! + Set smearTheseSensors = new HashSet(); + // I may eventually want this, so keep it here + private List _sensorsToSmear = new ArrayList(); + //List of Sensors + private List sensors = null; + private static final String SUBDETECTOR_NAME = "Tracker"; + private static Pattern layerPattern = Pattern.compile("L(\\d+)(t|b)"); + private static Pattern allPattern = Pattern.compile("L(\\all)"); + private Map fittedRawTrackerHitMap = new HashMap(); + private String fittedHitsCollectionName = "SVTFittedRawTrackerHits"; + private boolean _debug = false; + + Random r=new Random(); + + /// + + + + public void setSmearTimeFile(String smearFile){ + System.out.println("Setting SVT sensor time smearing file = "+smearFile); + this.smearTimeFile=smearFile; + } + + public void setDebug(boolean debug) { + this._debug = debug; + } + + public RawHitTimeSmearer() { + } + + @Override + public void detectorChanged(Detector detector) { + // Get the HpsSiSensor objects from the tracker detector element + sensors = detector.getSubdetector(SUBDETECTOR_NAME) + .getDetectorElement().findDescendants(HpsSiSensor.class); + // If the detector element had no sensors associated with it, throw + // an exception + if (sensors.size() == 0) + throw new RuntimeException(this.getClass().getName()+":: No sensors were found in this detector."); + if(_debug) + System.out.println(this.getClass().getName()+":: Reading in time smearing file = "+this.smearTimeFile); + Map mapOfSmearingTime=readSmearingFile(this.smearTimeFile); + if (mapOfSmearingTime.size()==0) + throw new RuntimeException(this.getClass().getName()+":: No sensors to smear???"); + + for (HpsSiSensor sensor: sensors){ + double smearingTime=-666.; + if(mapOfSmearingTime.containsKey(sensor.getName())){ + // found a sensor to smear...set up the object + if(_debug) + System.out.println(this.getClass().getName()+":: adding "+sensor.getName()+" with sigma time = "+mapOfSmearingTime.get(sensor.getName())); + smearingTime=mapOfSmearingTime.get(sensor.getName()); + } + if(smearingTime>0) + _sensorsToSmear.add(new SensorToSmear(sensor, smearingTime)); + } + if(_debug) + System.out.println(this.getClass().getName()+":: will smear cluster hits on "+_sensorsToSmear.size()+" sensors"); + + } + + @Override + public void process(EventHeader event) { + // get the sensor object of the cluster + for (SensorToSmear sensorToSmear : _sensorsToSmear){ + if(sensorToSmear.getSmearTimeSigma()>0){ + HpsSiSensor sensor=(HpsSiSensor)sensorToSmear.getSensor(); + List< LCRelation > fittedHits = sensor.getReadout().getHits(LCRelation.class); + for (LCRelation fittedHit : fittedHits) { + RawTrackerHit rth=FittedRawTrackerHit.getRawTrackerHit(fittedHit); + double oldTime=FittedRawTrackerHit.getT0(fittedHit); + double smearAmount=sensorToSmear.getRandomTimeSmear(); + double newTime=oldTime+smearAmount; + ((FittedRawTrackerHit)fittedHit).getShapeFitParameters().setT0(newTime); + if (_debug) + System.out.println("Smearing this hit Layer = " + sensor.getName() + +" smearTimeSigma= " + sensorToSmear.getSmearTimeSigma() + + " fitted hit time = " + FittedRawTrackerHit.getT0(fittedHit) + + " old time = " + oldTime + + " new time = " + newTime + + " amount smeared = " + smearAmount); + } + } + } + } + + @Override + public void endOfData() { + } + + public Map readSmearingFile(String smearFile){ + Map sensorNameSmearingMap = new HashMap(); + String infile = "/org/hps/recon/tracking/svtTimeAndPositionSmearing/" +smearFile; + if (_debug)System.out.println(this.getClass().getName()+"::Reading sensor smearing file " + infile); + InputStream inSensors = this.getClass().getResourceAsStream(infile); + BufferedReader reader = new BufferedReader(new InputStreamReader(inSensors)); + String line; + String delims = "\\s+";// this will split strings between one or more spaces + try { + while ((line = reader.readLine()) != null) { + String[] tokens = line.split(delims); + if (_debug) System.out.println("sensor name = " + tokens[0] + "; smearing = " + tokens[1]+" mm"); + String sensorName=tokens[0]; + Double smearingSigma=Double.parseDouble(tokens[1]); + sensorNameSmearingMap.put(sensorName,smearingSigma); + } + } catch (IOException ex) { + Logger.getLogger(this.getClass().getName()).log(Level.SEVERE, null, ex); + } + return sensorNameSmearingMap; + + } + + public class SensorToSmear { + + int _layer = 1; + HpsSiSensor _sensor = null; + double _smearTimeSigma = -666.0; // units of this is ns + + public SensorToSmear(HpsSiSensor sensor, double smearTime) { + _smearTimeSigma=smearTime; + _sensor = sensor; + } + + void setSmearTimeSigma(double smear) { + this._smearTimeSigma=smear; + } + + + double getSmearTimeSigma() { + return this._smearTimeSigma; + } + + HpsSiSensor getSensor(){ + return _sensor; + } + + boolean matchSensor(HpsSiSensor sensor) { + return _sensor == sensor; + } + + double getRandomTimeSmear(){ + if(this._smearTimeSigma>0) + return r.nextGaussian()*this._smearTimeSigma; + else + return 0.0; + } + } + + //Converts double array into Hep3Vector + private Hep3Vector toHep3(double[] arr) { + return new BasicHep3Vector(arr[0], arr[1], arr[2]); + } + + //Converts position into sensor frame + private Hep3Vector globalToSensor(Hep3Vector hitpos, HpsSiSensor sensor) { + SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + if (electrodes == null) { + electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON); + System.out.println("Charge Carrier is NULL"); + } + return electrodes.getGlobalToLocal().transformed(hitpos); + } + + //Converts position from local (sensor) to global frame + private Hep3Vector sensorToGlobal(Hep3Vector hitpos, HpsSiSensor sensor) { + SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + if (electrodes == null) { + electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON); + System.out.println("Charge Carrier is NULL"); + } + return electrodes.getLocalToGlobal().transformed(hitpos); + } + +} diff --git a/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java index 7bffec4a3..965990a0e 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java +++ b/tracking/src/main/java/org/hps/recon/tracking/mctweaking/StripHitSmearer.java @@ -13,6 +13,7 @@ import java.util.List; import java.util.Map; import java.util.Set; +import java.util.Random; import java.util.logging.Level; import java.util.logging.Logger; import java.util.regex.Matcher; @@ -23,94 +24,63 @@ import org.lcsim.event.EventHeader; import org.lcsim.event.RawTrackerHit; import org.lcsim.event.TrackerHit; +import org.lcsim.event.base.BaseTrackerHit; import org.lcsim.geometry.Detector; import org.lcsim.lcio.LCIOUtil; import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D; import org.lcsim.util.Driver; +import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType; +import hep.physics.matrix.SymmetricMatrix; /** * - * @author mgraham created 5/14/19 - * This driver will remove 1d strip clusters from the - * "StripClusterer_SiTrackerHitStrip1D" (default) - * collection based on a channel-based (data/mc) ratio file. - * - * Careful..the names of the ratio files are important! Format - * should be: - * _LX_top/bottom_stereo/axial_slot/hole.txt + * @author mgraham created 9/20/24 + * This driver will randomize the position of strip clusters + * within a specified gaussian sigma. * - * mg...this only works for L1 and L6 smearing at the moment - * ...unfortunately some of these things are hard coded + * Strip collection: + * "StripClusterer_SiTrackerHitStrip1D" (default) + * Important Settables: + * smear: the sigma of smearing gaussian + * + * mg...right now, I smear all hits in all layers by the same amount + * independent of nHits. */ public class StripHitSmearer extends Driver { + + String smearPositionFile="foobar"; //IMPORTANT...the layer, top/bottom/stereo/axial/slot/hole are derived from these names!!! - Set ratioFiles = new HashSet(); + Set smearTheseSensors = new HashSet(); + // I may eventually want this, so keep it here private List _sensorsToSmear = new ArrayList(); //List of Sensors private List sensors = null; private static final String SUBDETECTOR_NAME = "Tracker"; private static Pattern layerPattern = Pattern.compile("L(\\d+)(t|b)"); + private static Pattern allPattern = Pattern.compile("L(\\all)"); private String stripHitInputCollectionName = "StripClusterer_SiTrackerHitStrip1D"; private boolean _debug = false; - // instead of just removing hit at a given channel, remove all hits within Nsigma - private boolean removeHitsWithinNSig = false; - // used sigma-weighted ratios - private boolean useSigmaWeightedRatios = false; - double nSig = 5; - double sigmaUL1 = 0.1; //100 microns...this is roughly the 5-hit track projection error in U for layer 1, plotted in SVTHitLevelPlots - double sigmaUL6 = 0.25; //250 microns...this is roughly the 5-hit track projection error in U for layer 6*1.5 because pull is too wide, plotted in SVTHitLevelPlots - // double sigmaUL6 = 0.05; //50 microns...this is narrow just to check the effect. - // double firstSensorSmearFactor=3.0; - /// double secondSensorSmearFactor=2.0; - double firstSensorSmearFactor = 1.0; - double secondSensorSmearFactor = 1.0; - //these are just used for debugging... - int checkHitsChannel = 4; - int checkHitsTotal = 0; - int checkHitsPassed = 0; - double checkHitsRatio = 0; - double maxLayer = 12; + Random r=new Random(); + /// private List siClusters = new ArrayList(); - private Map _siClustersAcceptMap = new HashMap(); - private Map _finalSiClustersAcceptMap = new HashMap(); + // private Map _siClustersAcceptMap = new HashMap(); + //private Map _finalSiClustersAcceptMap = new HashMap(); - private Map _smearingFractionsL1 = new HashMap(); - private Map _smearingFractionsL6 = new HashMap(); - public void setRatioFiles(String[] ratioNames) { - System.out.println("Setting ratio files!!! " + ratioNames[0]); - this.ratioFiles = new HashSet(Arrays.asList(ratioNames)); + public void setSmearPositionFile(String smearFile){ + System.out.println("Setting SVT sensor position smearing file = "+smearFile); + this.smearPositionFile=smearFile; } - + public void setDebug(boolean debug) { this._debug = debug; } - public void setRemoveHitsWithinNSig(boolean removeHitsWithinNSig) { - this.removeHitsWithinNSig = removeHitsWithinNSig; - } - - public void setUseSigmaWeightedRatios(boolean useSigmaWeightedRatios) { - this.useSigmaWeightedRatios = useSigmaWeightedRatios; - } - - public void setL1Sigma(double l1sigma) { - this.sigmaUL1 = l1sigma; - } - - public void setL6Sigma(double l6sigma) { - this.sigmaUL1 = l6sigma; - } - - public void setNSigma(double nSig) { - this.nSig = nSig; - } - public StripHitSmearer() { } @@ -122,50 +92,33 @@ public void detectorChanged(Detector detector) { // If the detector element had no sensors associated with it, throw // an exception if (sensors.size() == 0) - throw new RuntimeException("No sensors were found in this detector."); - - //parse the ratio names and register sensors to kill - String delims = "_|\\.";// this will split strings between one "_" or "." - for (String ratioFile : ratioFiles) { - System.out.println("StripHitSmearer::Using this ratioFile: " + ratioFile); - int layer = -1; - boolean top = false; - boolean stereo = false; - boolean slot = false; - System.out.println("Parsing ratioFile = " + ratioFile); - String[] tokens = ratioFile.split(delims); - Matcher m = layerPattern.matcher(tokens[1]); - if (m.find()) { - layer = Integer.parseInt(m.group(1)); - if (m.group(2).matches("t")) - top = true; - } else { - System.out.println("Couldn't find layer number!!! " + ratioFile); - continue; - } -// if (tokens[2].matches("top")) -// top = true; -// System.out.println(tokens[2]+" "+tokens[3]+" "+tokens[4]); - if (tokens[2].matches("stereo")) - stereo = true; - if (tokens[3].matches("slot")) - slot = true; - System.out.println("StripHitSmearer::Smearing this: " - + "layer = " + layer + "; top = " + top + "; stereo = " + stereo - + "; slot = " + slot); - this.registerSensor(layer, top, stereo, slot, ratioFile); - } - System.out.println("filling smearing fractions"); - _smearingFractionsL1 = this.fillSmearedFractions(sigmaUL1, 1); - _smearingFractionsL6 = this.fillSmearedFractions(sigmaUL6, 12); + throw new RuntimeException(this.getClass().getName()+":: No sensors were found in this detector."); + if(_debug) + System.out.println(this.getClass().getName()+":: Reading in position smearing file = "+this.smearPositionFile); + Map mapOfSmearingPosition=readSmearingFile(this.smearPositionFile); + if (mapOfSmearingPosition.size()==0) + throw new RuntimeException(this.getClass().getName()+":: No sensors to smear???"); + + for (HpsSiSensor sensor: sensors){ + double smearingPosition=-666.; + if(mapOfSmearingPosition.containsKey(sensor.getName())){ + // found a sensor to smear...set up the object + if(_debug) + System.out.println(this.getClass().getName()+":: adding "+sensor.getName()+" with sigma = "+mapOfSmearingPosition.get(sensor.getName())); + smearingPosition=mapOfSmearingPosition.get(sensor.getName()); + // _sensorsToSmear.add(new SensorToSmear(sensor,mapOfSmearingPosition.get(sensor.getName()))); + } + if(smearingPosition>0) + _sensorsToSmear.add(new SensorToSmear(sensor,smearingPosition)); + } + if(_debug) + System.out.println(this.getClass().getName()+":: will smear cluster hits on "+_sensorsToSmear.size()+" sensors"); + } @Override public void process(EventHeader event) { - // System.out.println("In process of SVTHitSmearer"); - - int nhitsremoved = 0; - _siClustersAcceptMap.clear(); + List smearedSiClusters=new ArrayList(); if (event.hasItem(stripHitInputCollectionName)) siClusters = (List) event.get(stripHitInputCollectionName); else { @@ -175,371 +128,162 @@ public void process(EventHeader event) { if (_debug) System.out.println("Number of SiClusters Found = " + siClusters.size()); int oldClusterListSize = siClusters.size(); - for (TrackerHit siCluster : siClusters) { - boolean passHit = true; - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) siCluster.getRawHits().get(0)).getDetectorElement(); + for (TrackerHit siCluster : siClusters) { // Looping over all clusters + boolean smearedHit = false; + // get the sensor object of the cluster + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) siCluster.getRawHits().get(0)).getDetectorElement(); + Hep3Vector newPos=toHep3(siCluster.getPosition());//set the "new" variables to the original cluster position for (SensorToSmear sensorToSmear : _sensorsToSmear) - if (sensorToSmear.matchSensor(sensor)) { - //ok, get hit channel and kill or not - Hep3Vector pos = globalToSensor(toHep3(siCluster.getPosition()), sensor); - int chan = getChan(pos, sensor); - - double ratio = 0; - if (useSigmaWeightedRatios) - ratio = this.getSmearedRatio(chan, sensorToSmear); - else - ratio = sensorToSmear.getRatio(chan); - double killFactor = (1-ratio)*sensorToSmear.getScaleSmearFactor(); - if (_debug) - System.out.println("Found a hit on a sensor to kill!!! " + sensor.getName() + " Layer = " + sensor.getLayerNumber() - + " isTop? " + sensorToSmear.getIsTop() + " isStereo? " + sensorToSmear.getIsStereo() - + " isSlot? " + sensorToSmear.getIsSlot() + " channel = " + chan + " ratio = " + ratio); - - ratio = 1-killFactor; - if (ratio != -666) { - double random = Math.random(); //throw a random number to see if this hit should be rejected - if (random > ratio) { - passHit = false; - if (_debug) - System.out.println("Smearing this hit Layer = " + sensor.getLayerNumber() - + " isTop? " + sensorToSmear.getIsTop() + "; isStereo? " + sensorToSmear.getIsStereo() - + " isSlot? " + sensorToSmear.getIsSlot() +" scaleSmearFactor= " + sensorToSmear.getScaleSmearFactor() - + " channel = " + chan + " kill ratio = " + ratio + " random = " + random); - nhitsremoved++; - } - } - if (chan == checkHitsChannel) { - checkHitsTotal++; - if (passHit) - checkHitsPassed++; - checkHitsRatio = ratio; - } - - } - _siClustersAcceptMap.put(siCluster, passHit); + //get the sensorToSmear object in order to look up smearing size + if (sensorToSmear.matchSensor(sensor)) { + if(sensorToSmear.getSmearPositionSigma()>0){ + //get the local u position + Hep3Vector pos = globalToSensor(toHep3(siCluster.getPosition()), sensor); + + //this gives the amount to move the cluster based on the sigma assigned to the sensor + double smearAmount=sensorToSmear.getRandomPositionSmear(); + double uPosNew=pos.x()+smearAmount; + //get the new global position + newPos=sensorToGlobal(toHep3(new double[]{uPosNew,pos.y(),pos.z()}),sensor); + //make a copy of cluster and set the new position + // SiTrackerHitStrip1D smearedTrackerHit=new SiTrackerHitStrip1D(newPos, new SymmetricMatrix(3,siCluster.getCovMatrix(),true), + //siCluster.getdEdx(),siCluster.getTime(),(List)siCluster.getRawHits(),TrackerHitType.decoded(siCluster.getType())); + + if (_debug) + System.out.println("Smearing this hit Layer = " + sensor.getName() + +" smearPositionSigma= " + sensorToSmear.getSmearPositionSigma() + + " old u = " + pos.x() + + " new u = " + uPosNew + + " amount smeared = " + smearAmount); + if(_debug){ + System.out.println("original position = "+toHep3(siCluster.getPosition()).toString()); + System.out.println("global og position = "+pos.toString()); + } + // set the new position + // smearedTrackerHit.setPosition(new double[]{newPos.x(),newPos.y(),newPos.z()}); + + // add the smeared hit to list and mark it + // smearedSiClusters.add(smearedTrackerHit); + smearedHit=true; + } + } + //make a smeared hit...even if the sensor hit was not smeared, use newPos/newTime since they were filled with original pos/time above + SiTrackerHitStrip1D smearedTrackerHit= + new SiTrackerHitStrip1D(newPos, new SymmetricMatrix(3,siCluster.getCovMatrix(),true), + siCluster.getdEdx(),siCluster.getTime(),(List)siCluster.getRawHits(), + TrackerHitType.decoded(siCluster.getType())); + + smearedSiClusters.add(smearedTrackerHit); + } -// if (_debug) - List tmpClusterList = getFinalHits(_siClustersAcceptMap); - // if (_debug) - if (_debug)System.out.println("New Cluster List Has " + tmpClusterList.size() + "; old List had " + oldClusterListSize); - if (_debug)System.out.println("number of hits removed for this event = " + nhitsremoved); + if (_debug)System.out.println("New Cluster List Has " + smearedSiClusters.size() + "; old List had " + oldClusterListSize); // TODO flag not used int flag = LCIOUtil.bitSet(0, 31, true); // Turn on 64-bit cell ID. event.remove(this.stripHitInputCollectionName); - event.put(this.stripHitInputCollectionName, tmpClusterList, SiTrackerHitStrip1D.class, 0, toString()); - + event.put(this.stripHitInputCollectionName, smearedSiClusters, SiTrackerHitStrip1D.class, 0, toString()); } @Override public void endOfData() { - System.out.println("StripHitSmearer::endOfData channel = " + checkHitsChannel - + " ratio = " + checkHitsRatio + "; hits pass = " + checkHitsPassed - + " hits total = " + checkHitsTotal + " passRatio = " + ((double) checkHitsPassed) / checkHitsTotal); - } - - public void registerSensor(int layer, boolean isTop, boolean isStereo, boolean isSlot, String ratioFile) { - if (_debug) { - System.out.println("================ Registering New SensorToSmear =============== "); - System.out.println(" " - + "layer = " + layer + "; top = " + isTop + "; stereo = " + isStereo - + "; slot = " + isSlot + "; ratio file = " + ratioFile); - } - SensorToSmear newSensor = new SensorToSmear(layer, isTop, isStereo, isSlot, ratioFile); - System.out.println("newSensor isTop " + newSensor.getIsTop()); - if (newSensor.getIsTop() && !(newSensor.getIsStereo()) && newSensor.getLayer()==1) - newSensor.setScaleSmearFactor(firstSensorSmearFactor) ; - if (newSensor.getIsTop() && newSensor.getIsStereo() && newSensor.getLayer()==1) - newSensor.setScaleSmearFactor(secondSensorSmearFactor); - - if (!(newSensor.getIsTop()) && newSensor.getIsStereo() && newSensor.getLayer()==1) - newSensor.setScaleSmearFactor(firstSensorSmearFactor); - if (!(newSensor.getIsTop()) && !(newSensor.getIsStereo()) && newSensor.getLayer()==1) - newSensor.setScaleSmearFactor(secondSensorSmearFactor); - if (_debug) { - System.out.println(" channel 10 ratio for newSensor = "+newSensor.getRatio(10)); - System.out.println("================ Done With New SensorToSmear =============== "); - } - _sensorsToSmear.add(newSensor); - } - - //Return the HpsSiSensor for a given top/bottom track, layer, axial/stereo, and slot/hole - private HpsSiSensor getSensor(int layer, boolean isTop, boolean isAxial, boolean isHole) { - for (HpsSiSensor sensor : sensors) { - int senselayer = (sensor.getLayerNumber() + 1) / 2; - if (senselayer != layer) - continue; - if ((isTop && !sensor.isTopLayer()) || (!isTop && sensor.isTopLayer())) - continue; - if ((isAxial && !sensor.isAxial()) || (!isAxial && sensor.isAxial())) - continue; - if (layer < 4 && layer > 0) {//only for pre-2019 runs...fix this... - if (_debug)System.out.println("Matching with sensor = " + sensor.getName()); - return sensor; - } - else { - if ((!sensor.getSide().matches("ELECTRON") && isHole) || (sensor.getSide().matches("ELECTRON") && !isHole)) - continue; - - if (_debug)System.out.println("Matching with sensor = " + sensor.getName()); - return sensor; - } - } - return null; - } - - private List getFinalHits(Map _initialSiClustersAcceptMap) { - _finalSiClustersAcceptMap.clear(); - List tmpClusterList = new ArrayList(); - for (TrackerHit hit : _initialSiClustersAcceptMap.keySet()) { - boolean isAccept = _initialSiClustersAcceptMap.get(hit); - if (removeHitsWithinNSig && !isAccept) { // this hit got thrown out; lets throw out all hits within N-sigma - if (_debug) - System.out.println("Found a killed hit...now check for hits around it"); - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) hit.getRawHits().get(0)).getDetectorElement(); - Hep3Vector pos = globalToSensor(toHep3(hit.getPosition()), sensor); - int chan = getChan(pos, sensor); - double sigmaU = sigmaUL1; - if (sensor.getLayerNumber() > 10) //remember, layer number is sensor number (1-12) - sigmaU = sigmaUL6; - if (_debug) - System.out.println("Layer number = " + sensor.getLayerNumber() + " sigma = " + sigmaU); - Hep3Vector minPos = new BasicHep3Vector(pos.x() - nSig * sigmaU, pos.y(), pos.z()); - Hep3Vector maxPos = new BasicHep3Vector(pos.x() + nSig * sigmaU, pos.y(), pos.z()); - int minChan = getChan(minPos, sensor); - int maxChan = getChan(maxPos, sensor); - if (minChan > maxChan) { - int tmp = minChan; - minChan = maxChan; - maxChan = tmp; - } - if (_debug) - System.out.println("hit channel = " + chan + "; minChan = " + minChan + "; maxChan = " + maxChan); - //loop through initial clusters again - for (TrackerHit testhit : _initialSiClustersAcceptMap.keySet()) { - HpsSiSensor testsensor = (HpsSiSensor) ((RawTrackerHit) testhit.getRawHits().get(0)).getDetectorElement(); - if (sensor != testsensor) - continue; - Hep3Vector testpos = globalToSensor(toHep3(testhit.getPosition()), sensor); - int testchan = getChan(testpos, sensor); - if (_debug) - System.out.println("test channel = " + testchan); - if (testchan > maxChan || testchan < minChan) - continue; - //in range, set accept to false - if (_debug) - System.out.println("Found a hit within range of killed hit...killing it too!"); - _finalSiClustersAcceptMap.put(testhit, false); - } - _finalSiClustersAcceptMap.put(hit, isAccept); //add the original cluster to map (shouldn't matter since isAccept==false - } else if (!_finalSiClustersAcceptMap.containsKey(hit)) - _finalSiClustersAcceptMap.put(hit, isAccept); - } - //fill tmpClusterList with accepted hits - for (TrackerHit hit : _finalSiClustersAcceptMap.keySet()) - if (_finalSiClustersAcceptMap.get(hit)) - tmpClusterList.add(hit); - return tmpClusterList; - } - - private double getSmearedRatio(int seedChan, SensorToSmear sensorToSmear) { - double smearedRatio = 0; - double integralSum = 0; - Map _smearingFractions = new HashMap(); - if (sensorToSmear.getLayer() == 1) - _smearingFractions = _smearingFractionsL1; - else if (sensorToSmear.getLayer() == 6) - _smearingFractions = _smearingFractionsL6; - for (Map.Entry entry : _smearingFractions.entrySet()) { - int chan = seedChan + entry.getKey(); - if (chan < 0 || chan > sensorToSmear._sensor.getNumberOfChannels() - 1) - continue; - double chanIntegral = entry.getValue(); - double ratio = sensorToSmear.getRatio(chan); - if (ratio > 1.0 || ratio == -666) - ratio = 1.0; -// System.out.println("channel = " + chan + " ratio = " - // + ratio + "; smearedRatio = " + ratio * chanIntegral); - smearedRatio += ratio * chanIntegral; - integralSum += chanIntegral; - } - if (integralSum == 0) - smearedRatio = 1.0; - else - smearedRatio /= integralSum; - if (_debug) - System.out.println("getSmearedRatio: seed ratio = " - + sensorToSmear.getRatio(seedChan) + "; smearedRatio = " + smearedRatio); - return smearedRatio; - } - - private Map fillSmearedFractions(double sigmaU, int layer) { - Map _smearingFractions = new HashMap(); - double pitch = 0; - for (HpsSiSensor s : sensors) - if (s.getLayerNumber() == layer) - pitch = s.getReadoutStripPitch(); - if (pitch == 0) { - System.out.println("Couldn't find layer = " + layer); - return _smearingFractions; - } - - double normPitch = sigmaU / pitch; - int nStrips = (int) Math.floor(this.nSig * normPitch); - double totalInt = 0; - for (int i = -nStrips; i < 1; i++) {//just do negative bins (and 0)..positive side is same - double minIntLimit = (i - 0.5) / normPitch; - double maxIntLimit = (i + 0.5) / normPitch; - double minInt = this.computeGaussInt(minIntLimit, 1000); - double maxInt = this.computeGaussInt(maxIntLimit, 1000); - totalInt += maxInt - minInt; - _smearingFractions.put(i, maxInt - minInt); - System.out.println("_smearingFractions chan = " + i + " fraction = " + _smearingFractions.get(i)); - - } - for (int i = 1; i < nStrips + 1; i++) {//positive side is same - totalInt += _smearingFractions.get(-1 * i); - _smearingFractions.put(i, _smearingFractions.get(-1 * i)); - System.out.println("_smearingFractions chan = " + i + " fraction = " + _smearingFractions.get(i)); - - } - System.out.println("Layer = " + layer + "; totalInt = " + totalInt); - return _smearingFractions; } + public Map readSmearingFile(String smearFile){ + Map sensorNameSmearingMap = new HashMap(); + String infile = "/org/hps/recon/tracking/svtTimeAndPositionSmearing/" +smearFile; + if (_debug)System.out.println(this.getClass().getName()+"::Reading sensor smearing file " + infile); + InputStream inSensors = this.getClass().getResourceAsStream(infile); + BufferedReader reader = new BufferedReader(new InputStreamReader(inSensors)); + String line; + String delims = "\\s+";// this will split strings between one or more spaces + try { + while ((line = reader.readLine()) != null) { + String[] tokens = line.split(delims); + if (_debug) System.out.println("sensor name = " + tokens[0] + "; smearing = " + tokens[1]+" mm"); + String sensorName=tokens[0]; + Double smearingSigma=Double.parseDouble(tokens[1]); + sensorNameSmearingMap.put(sensorName,smearingSigma); + } + } catch (IOException ex) { + Logger.getLogger(this.getClass().getName()).log(Level.SEVERE, null, ex); + } + return sensorNameSmearingMap; + + } + public class SensorToSmear { int _layer = 1; - boolean _isStereo = false; - boolean _isSlot = false; - boolean _isTop = false; - String _ratioFile = "foobarTopL1Stereo.txt"; HpsSiSensor _sensor = null; - double _scaleSmearFactor = 1.0; - Map _channelToRatioMap = new HashMap(); + double _smearPositionSigma = -666.0; // units of this is mm + double _smearTimeSigma = -666.0; // units of this is ns - public SensorToSmear(int layer, boolean isTop, boolean isStereo, boolean isSlot, String ratioFile) { - if (_debug)System.out.println("Making new SensorToSmear layer = " + layer); - _layer = layer; - _isTop = isTop; - _isStereo = isStereo; - _isSlot = isSlot; - _ratioFile = ratioFile; - _sensor = getSensor(layer, isTop, !isStereo, !isSlot); - readRatioFile(); + public SensorToSmear(HpsSiSensor sensor, double smearPos,double smearTime) { + _smearPositionSigma=smearPos; + _smearTimeSigma=smearTime; + _sensor = sensor; } - - void setScaleSmearFactor(double scale) { - this._scaleSmearFactor=scale; - } - - int getLayer() { - return _layer; - } - - boolean getIsStereo() { - return _isStereo; + public SensorToSmear(HpsSiSensor sensor, double smearPos) { + _smearPositionSigma=smearPos; + _sensor = sensor; } - - boolean getIsSlot() { - return _isSlot; + void setSmearPositionSigma(double smear) { + this._smearPositionSigma=smear; } - - boolean getIsTop() { - return _isTop; + + void setSmearTimeSigma(double smear) { + this._smearTimeSigma=smear; } - - String getRatioFile() { - return _ratioFile; + + double getSmearPositionSigma() { + return this._smearPositionSigma; } - - double getScaleSmearFactor() { - return(this._scaleSmearFactor); + + double getSmearTimeSigma() { + return this._smearTimeSigma; } - + boolean matchSensor(HpsSiSensor sensor) { return _sensor == sensor; } - private void readRatioFile() { - String infile = "/org/hps/recon/tracking/efficiencyCorrections/" + _ratioFile; - InputStream inRatios = this.getClass().getResourceAsStream(infile); - if (_debug)System.out.println("StripHitSmearer::Reading ratio file " + infile); - BufferedReader reader = new BufferedReader(new InputStreamReader(inRatios)); - String line; - String delims = "[ ]+";// this will split strings between one or more spaces - try { - while ((line = reader.readLine()) != null) { - String[] tokens = line.split(delims); - if (_debug) System.out.println("channel number = " + tokens[0] + "; ratio = " + tokens[1]); - _channelToRatioMap.put(Integer.parseInt(tokens[0]), Double.parseDouble(tokens[1])); - } - } catch (IOException ex) { - Logger.getLogger(StripHitSmearer.class.getName()).log(Level.SEVERE, null, ex); - } - } - - public double getRatio(int channel) { - if (_debug)System.out.println("SensorToSmear:: getRatio: layer = " + this.getLayer() + "; channel " + channel); - if (_channelToRatioMap.get(channel) == null) - return -666; - return _channelToRatioMap.get(channel); - } - + double getRandomPositionSmear(){ + if(this._smearPositionSigma>0) + return r.nextGaussian()*this._smearPositionSigma; + else + return 0.0; + } + double getRandomTimeSmear(){ + if(this._smearTimeSigma>0) + return r.nextGaussian()*this._smearTimeSigma; + else + return 0.0; + } } //Converts double array into Hep3Vector private Hep3Vector toHep3(double[] arr) { return new BasicHep3Vector(arr[0], arr[1], arr[2]); } - //Returns channel number of a given position in the sensor frame - /* - private int getChan(Hep3Vector pos, HpsSiSensor sensor) { - double readoutPitch = sensor.getReadoutStripPitch(); - int nChan = sensor.getNumberOfChannels(); - double height = readoutPitch * nChan; - return (int) ((height / 2 - pos.x()) / readoutPitch); - } - */ - - private int getChan(Hep3Vector pos, HpsSiSensor sensor) { - SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); - if (maxLayer>12 && sensor.getLayerNumber()<5) {//thin sensors - int row = electrodes.getRowNumber(pos); - int col = electrodes.getColumnNumber(pos); - return electrodes.getCellID(row,col); - }else{ - return electrodes.getCellID(pos); - } - } + //Converts position into sensor frame - private Hep3Vector globalToSensor(Hep3Vector trkpos, HpsSiSensor sensor) { + private Hep3Vector globalToSensor(Hep3Vector hitpos, HpsSiSensor sensor) { SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); if (electrodes == null) { electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON); System.out.println("Charge Carrier is NULL"); } - return electrodes.getGlobalToLocal().transformed(trkpos); + return electrodes.getGlobalToLocal().transformed(hitpos); } - //Computes gaussian integral numerically from -inf to nSig - // this should be replaced by a simple error function but I can't find one - // in normal java libraries!!! - private double computeGaussInt(double nSig, int nSteps) { - double mean = 0; - double sigma = 1; - double dx = sigma * nSig / (double) nSteps; - double integral = 0; - for (int i = 0; i < nSteps; i++) { - double x = dx * (i + 0.5) + mean; - integral += dx * Gauss(x, mean, sigma); + //Converts position from local (sensor) to global frame + private Hep3Vector sensorToGlobal(Hep3Vector hitpos, HpsSiSensor sensor) { + SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + if (electrodes == null) { + electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON); + System.out.println("Charge Carrier is NULL"); } - return integral + 0.5; - } - - //Gaussian function - private double Gauss(double x, double mean, double sigma) { - return 1 / (Math.sqrt(2 * Math.PI * Math.pow(sigma, 2))) * Math.exp(-Math.pow(x - mean, 2) / (2 * Math.pow(sigma, 2))); + return electrodes.getLocalToGlobal().transformed(hitpos); } } diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt008mm.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt008mm.txt new file mode 100644 index 000000000..a1f95db04 --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt008mm.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 0.008 +module_L1b_halfmodule_stereo_sensor0 0.008 +module_L2b_halfmodule_axial_sensor0 0.008 +module_L2b_halfmodule_stereo_sensor0 0.008 +module_L3b_halfmodule_axial_sensor0 0.008 +module_L3b_halfmodule_stereo_sensor0 0.008 +module_L4b_halfmodule_axial_hole_sensor0 0.008 +module_L4b_halfmodule_axial_slot_sensor0 0.008 +module_L4b_halfmodule_stereo_hole_sensor0 0.008 +module_L4b_halfmodule_stereo_slot_sensor0 0.008 +module_L5b_halfmodule_axial_hole_sensor0 0.008 +module_L5b_halfmodule_axial_slot_sensor0 0.008 +module_L5b_halfmodule_stereo_hole_sensor0 0.008 +module_L5b_halfmodule_stereo_slot_sensor0 0.008 +module_L6b_halfmodule_axial_hole_sensor0 0.008 +module_L6b_halfmodule_axial_slot_sensor0 0.008 +module_L6b_halfmodule_stereo_hole_sensor0 0.008 +module_L6b_halfmodule_stereo_slot_sensor0 0.008 +module_L1t_halfmodule_axial_sensor0 0.008 +module_L1t_halfmodule_stereo_sensor0 0.008 +module_L2t_halfmodule_axial_sensor0 0.008 +module_L2t_halfmodule_stereo_sensor0 0.008 +module_L3t_halfmodule_axial_sensor0 0.008 +module_L3t_halfmodule_stereo_sensor0 0.008 +module_L4t_halfmodule_axial_hole_sensor0 0.008 +module_L4t_halfmodule_axial_slot_sensor0 0.008 +module_L4t_halfmodule_stereo_hole_sensor0 0.008 +module_L4t_halfmodule_stereo_slot_sensor0 0.008 +module_L5t_halfmodule_axial_hole_sensor0 0.008 +module_L5t_halfmodule_axial_slot_sensor0 0.008 +module_L5t_halfmodule_stereo_hole_sensor0 0.008 +module_L5t_halfmodule_stereo_slot_sensor0 0.008 +module_L6t_halfmodule_axial_hole_sensor0 0.008 +module_L6t_halfmodule_axial_slot_sensor0 0.008 +module_L6t_halfmodule_stereo_hole_sensor0 0.008 +module_L6t_halfmodule_stereo_slot_sensor0 0.008 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt01mm.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt01mm.txt new file mode 100644 index 000000000..48a04d84b --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt01mm.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 0.01 +module_L1b_halfmodule_stereo_sensor0 0.01 +module_L2b_halfmodule_axial_sensor0 0.01 +module_L2b_halfmodule_stereo_sensor0 0.01 +module_L3b_halfmodule_axial_sensor0 0.01 +module_L3b_halfmodule_stereo_sensor0 0.01 +module_L4b_halfmodule_axial_hole_sensor0 0.01 +module_L4b_halfmodule_axial_slot_sensor0 0.01 +module_L4b_halfmodule_stereo_hole_sensor0 0.01 +module_L4b_halfmodule_stereo_slot_sensor0 0.01 +module_L5b_halfmodule_axial_hole_sensor0 0.01 +module_L5b_halfmodule_axial_slot_sensor0 0.01 +module_L5b_halfmodule_stereo_hole_sensor0 0.01 +module_L5b_halfmodule_stereo_slot_sensor0 0.01 +module_L6b_halfmodule_axial_hole_sensor0 0.01 +module_L6b_halfmodule_axial_slot_sensor0 0.01 +module_L6b_halfmodule_stereo_hole_sensor0 0.01 +module_L6b_halfmodule_stereo_slot_sensor0 0.01 +module_L1t_halfmodule_axial_sensor0 0.01 +module_L1t_halfmodule_stereo_sensor0 0.01 +module_L2t_halfmodule_axial_sensor0 0.01 +module_L2t_halfmodule_stereo_sensor0 0.01 +module_L3t_halfmodule_axial_sensor0 0.01 +module_L3t_halfmodule_stereo_sensor0 0.01 +module_L4t_halfmodule_axial_hole_sensor0 0.01 +module_L4t_halfmodule_axial_slot_sensor0 0.01 +module_L4t_halfmodule_stereo_hole_sensor0 0.01 +module_L4t_halfmodule_stereo_slot_sensor0 0.01 +module_L5t_halfmodule_axial_hole_sensor0 0.01 +module_L5t_halfmodule_axial_slot_sensor0 0.01 +module_L5t_halfmodule_stereo_hole_sensor0 0.01 +module_L5t_halfmodule_stereo_slot_sensor0 0.01 +module_L6t_halfmodule_axial_hole_sensor0 0.01 +module_L6t_halfmodule_axial_slot_sensor0 0.01 +module_L6t_halfmodule_stereo_hole_sensor0 0.01 +module_L6t_halfmodule_stereo_slot_sensor0 0.01 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt0mm.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt0mm.txt new file mode 100644 index 000000000..0ee5ea79a --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt0mm.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 0.0 +module_L1b_halfmodule_stereo_sensor0 0.0 +module_L2b_halfmodule_axial_sensor0 0.0 +module_L2b_halfmodule_stereo_sensor0 0.0 +module_L3b_halfmodule_axial_sensor0 0.0 +module_L3b_halfmodule_stereo_sensor0 0.0 +module_L4b_halfmodule_axial_hole_sensor0 0.0 +module_L4b_halfmodule_axial_slot_sensor0 0.0 +module_L4b_halfmodule_stereo_hole_sensor0 0.0 +module_L4b_halfmodule_stereo_slot_sensor0 0.0 +module_L5b_halfmodule_axial_hole_sensor0 0.0 +module_L5b_halfmodule_axial_slot_sensor0 0.0 +module_L5b_halfmodule_stereo_hole_sensor0 0.0 +module_L5b_halfmodule_stereo_slot_sensor0 0.0 +module_L6b_halfmodule_axial_hole_sensor0 0.0 +module_L6b_halfmodule_axial_slot_sensor0 0.0 +module_L6b_halfmodule_stereo_hole_sensor0 0.0 +module_L6b_halfmodule_stereo_slot_sensor0 0.0 +module_L1t_halfmodule_axial_sensor0 0.0 +module_L1t_halfmodule_stereo_sensor0 0.0 +module_L2t_halfmodule_axial_sensor0 0.0 +module_L2t_halfmodule_stereo_sensor0 0.0 +module_L3t_halfmodule_axial_sensor0 0.0 +module_L3t_halfmodule_stereo_sensor0 0.0 +module_L4t_halfmodule_axial_hole_sensor0 0.0 +module_L4t_halfmodule_axial_slot_sensor0 0.0 +module_L4t_halfmodule_stereo_hole_sensor0 0.0 +module_L4t_halfmodule_stereo_slot_sensor0 0.0 +module_L5t_halfmodule_axial_hole_sensor0 0.0 +module_L5t_halfmodule_axial_slot_sensor0 0.0 +module_L5t_halfmodule_stereo_hole_sensor0 0.0 +module_L5t_halfmodule_stereo_slot_sensor0 0.0 +module_L6t_halfmodule_axial_hole_sensor0 0.0 +module_L6t_halfmodule_axial_slot_sensor0 0.0 +module_L6t_halfmodule_stereo_hole_sensor0 0.0 +module_L6t_halfmodule_stereo_slot_sensor0 0.0 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt1mm.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt1mm.txt new file mode 100644 index 000000000..5c794e237 --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/exampleSmearing-0pt1mm.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 0.1 +module_L1b_halfmodule_stereo_sensor0 0.1 +module_L2b_halfmodule_axial_sensor0 0.1 +module_L2b_halfmodule_stereo_sensor0 0.1 +module_L3b_halfmodule_axial_sensor0 0.1 +module_L3b_halfmodule_stereo_sensor0 0.1 +module_L4b_halfmodule_axial_hole_sensor0 0.1 +module_L4b_halfmodule_axial_slot_sensor0 0.1 +module_L4b_halfmodule_stereo_hole_sensor0 0.1 +module_L4b_halfmodule_stereo_slot_sensor0 0.1 +module_L5b_halfmodule_axial_hole_sensor0 0.1 +module_L5b_halfmodule_axial_slot_sensor0 0.1 +module_L5b_halfmodule_stereo_hole_sensor0 0.1 +module_L5b_halfmodule_stereo_slot_sensor0 0.1 +module_L6b_halfmodule_axial_hole_sensor0 0.1 +module_L6b_halfmodule_axial_slot_sensor0 0.1 +module_L6b_halfmodule_stereo_hole_sensor0 0.1 +module_L6b_halfmodule_stereo_slot_sensor0 0.1 +module_L1t_halfmodule_axial_sensor0 0.1 +module_L1t_halfmodule_stereo_sensor0 0.1 +module_L2t_halfmodule_axial_sensor0 0.1 +module_L2t_halfmodule_stereo_sensor0 0.1 +module_L3t_halfmodule_axial_sensor0 0.1 +module_L3t_halfmodule_stereo_sensor0 0.1 +module_L4t_halfmodule_axial_hole_sensor0 0.1 +module_L4t_halfmodule_axial_slot_sensor0 0.1 +module_L4t_halfmodule_stereo_hole_sensor0 0.1 +module_L4t_halfmodule_stereo_slot_sensor0 0.1 +module_L5t_halfmodule_axial_hole_sensor0 0.1 +module_L5t_halfmodule_axial_slot_sensor0 0.1 +module_L5t_halfmodule_stereo_hole_sensor0 0.1 +module_L5t_halfmodule_stereo_slot_sensor0 0.1 +module_L6t_halfmodule_axial_hole_sensor0 0.1 +module_L6t_halfmodule_axial_slot_sensor0 0.1 +module_L6t_halfmodule_stereo_hole_sensor0 0.1 +module_L6t_halfmodule_stereo_slot_sensor0 0.1 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns-4nsL0.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns-4nsL0.txt new file mode 100644 index 000000000..2f82f048c --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns-4nsL0.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 4.0 +module_L1b_halfmodule_stereo_sensor0 4.0 +module_L2b_halfmodule_axial_sensor0 2.0 +module_L2b_halfmodule_stereo_sensor0 2.0 +module_L3b_halfmodule_axial_sensor0 2.0 +module_L3b_halfmodule_stereo_sensor0 2.0 +module_L4b_halfmodule_axial_hole_sensor0 2.0 +module_L4b_halfmodule_axial_slot_sensor0 2.0 +module_L4b_halfmodule_stereo_hole_sensor0 2.0 +module_L4b_halfmodule_stereo_slot_sensor0 2.0 +module_L5b_halfmodule_axial_hole_sensor0 2.0 +module_L5b_halfmodule_axial_slot_sensor0 2.0 +module_L5b_halfmodule_stereo_hole_sensor0 2.0 +module_L5b_halfmodule_stereo_slot_sensor0 2.0 +module_L6b_halfmodule_axial_hole_sensor0 2.0 +module_L6b_halfmodule_axial_slot_sensor0 2.0 +module_L6b_halfmodule_stereo_hole_sensor0 2.0 +module_L6b_halfmodule_stereo_slot_sensor0 2.0 +module_L1t_halfmodule_axial_sensor0 4.0 +module_L1t_halfmodule_stereo_sensor0 4.0 +module_L2t_halfmodule_axial_sensor0 2.0 +module_L2t_halfmodule_stereo_sensor0 2.0 +module_L3t_halfmodule_axial_sensor0 2.0 +module_L3t_halfmodule_stereo_sensor0 2.0 +module_L4t_halfmodule_axial_hole_sensor0 2.0 +module_L4t_halfmodule_axial_slot_sensor0 2.0 +module_L4t_halfmodule_stereo_hole_sensor0 2.0 +module_L4t_halfmodule_stereo_slot_sensor0 2.0 +module_L5t_halfmodule_axial_hole_sensor0 2.0 +module_L5t_halfmodule_axial_slot_sensor0 2.0 +module_L5t_halfmodule_stereo_hole_sensor0 2.0 +module_L5t_halfmodule_stereo_slot_sensor0 2.0 +module_L6t_halfmodule_axial_hole_sensor0 2.0 +module_L6t_halfmodule_axial_slot_sensor0 2.0 +module_L6t_halfmodule_stereo_hole_sensor0 2.0 +module_L6t_halfmodule_stereo_slot_sensor0 2.0 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns.txt new file mode 100644 index 000000000..836987547 --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-2ns.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 2.0 +module_L1b_halfmodule_stereo_sensor0 2.0 +module_L2b_halfmodule_axial_sensor0 2.0 +module_L2b_halfmodule_stereo_sensor0 2.0 +module_L3b_halfmodule_axial_sensor0 2.0 +module_L3b_halfmodule_stereo_sensor0 2.0 +module_L4b_halfmodule_axial_hole_sensor0 2.0 +module_L4b_halfmodule_axial_slot_sensor0 2.0 +module_L4b_halfmodule_stereo_hole_sensor0 2.0 +module_L4b_halfmodule_stereo_slot_sensor0 2.0 +module_L5b_halfmodule_axial_hole_sensor0 2.0 +module_L5b_halfmodule_axial_slot_sensor0 2.0 +module_L5b_halfmodule_stereo_hole_sensor0 2.0 +module_L5b_halfmodule_stereo_slot_sensor0 2.0 +module_L6b_halfmodule_axial_hole_sensor0 2.0 +module_L6b_halfmodule_axial_slot_sensor0 2.0 +module_L6b_halfmodule_stereo_hole_sensor0 2.0 +module_L6b_halfmodule_stereo_slot_sensor0 2.0 +module_L1t_halfmodule_axial_sensor0 2.0 +module_L1t_halfmodule_stereo_sensor0 2.0 +module_L2t_halfmodule_axial_sensor0 2.0 +module_L2t_halfmodule_stereo_sensor0 2.0 +module_L3t_halfmodule_axial_sensor0 2.0 +module_L3t_halfmodule_stereo_sensor0 2.0 +module_L4t_halfmodule_axial_hole_sensor0 2.0 +module_L4t_halfmodule_axial_slot_sensor0 2.0 +module_L4t_halfmodule_stereo_hole_sensor0 2.0 +module_L4t_halfmodule_stereo_slot_sensor0 2.0 +module_L5t_halfmodule_axial_hole_sensor0 2.0 +module_L5t_halfmodule_axial_slot_sensor0 2.0 +module_L5t_halfmodule_stereo_hole_sensor0 2.0 +module_L5t_halfmodule_stereo_slot_sensor0 2.0 +module_L6t_halfmodule_axial_hole_sensor0 2.0 +module_L6t_halfmodule_axial_slot_sensor0 2.0 +module_L6t_halfmodule_stereo_hole_sensor0 2.0 +module_L6t_halfmodule_stereo_slot_sensor0 2.0 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-5ns.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-5ns.txt new file mode 100644 index 000000000..3727b36bc --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-5ns.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 5.0 +module_L1b_halfmodule_stereo_sensor0 5.0 +module_L2b_halfmodule_axial_sensor0 5.0 +module_L2b_halfmodule_stereo_sensor0 5.0 +module_L3b_halfmodule_axial_sensor0 5.0 +module_L3b_halfmodule_stereo_sensor0 5.0 +module_L4b_halfmodule_axial_hole_sensor0 5.0 +module_L4b_halfmodule_axial_slot_sensor0 5.0 +module_L4b_halfmodule_stereo_hole_sensor0 5.0 +module_L4b_halfmodule_stereo_slot_sensor0 5.0 +module_L5b_halfmodule_axial_hole_sensor0 5.0 +module_L5b_halfmodule_axial_slot_sensor0 5.0 +module_L5b_halfmodule_stereo_hole_sensor0 5.0 +module_L5b_halfmodule_stereo_slot_sensor0 5.0 +module_L6b_halfmodule_axial_hole_sensor0 5.0 +module_L6b_halfmodule_axial_slot_sensor0 5.0 +module_L6b_halfmodule_stereo_hole_sensor0 5.0 +module_L6b_halfmodule_stereo_slot_sensor0 5.0 +module_L1t_halfmodule_axial_sensor0 5.0 +module_L1t_halfmodule_stereo_sensor0 5.0 +module_L2t_halfmodule_axial_sensor0 5.0 +module_L2t_halfmodule_stereo_sensor0 5.0 +module_L3t_halfmodule_axial_sensor0 5.0 +module_L3t_halfmodule_stereo_sensor0 5.0 +module_L4t_halfmodule_axial_hole_sensor0 5.0 +module_L4t_halfmodule_axial_slot_sensor0 5.0 +module_L4t_halfmodule_stereo_hole_sensor0 5.0 +module_L4t_halfmodule_stereo_slot_sensor0 5.0 +module_L5t_halfmodule_axial_hole_sensor0 5.0 +module_L5t_halfmodule_axial_slot_sensor0 5.0 +module_L5t_halfmodule_stereo_hole_sensor0 5.0 +module_L5t_halfmodule_stereo_slot_sensor0 5.0 +module_L6t_halfmodule_axial_hole_sensor0 5.0 +module_L6t_halfmodule_axial_slot_sensor0 5.0 +module_L6t_halfmodule_stereo_hole_sensor0 5.0 +module_L6t_halfmodule_stereo_slot_sensor0 5.0 diff --git a/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-almost-zero.txt b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-almost-zero.txt new file mode 100644 index 000000000..bfe6e9557 --- /dev/null +++ b/tracking/src/main/resources/org/hps/recon/tracking/svtTimeAndPositionSmearing/timeSmearing-almost-zero.txt @@ -0,0 +1,36 @@ +module_L1b_halfmodule_axial_sensor0 0.00001 +module_L1b_halfmodule_stereo_sensor0 0.00001 +module_L2b_halfmodule_axial_sensor0 0.00001 +module_L2b_halfmodule_stereo_sensor0 0.00001 +module_L3b_halfmodule_axial_sensor0 0.00001 +module_L3b_halfmodule_stereo_sensor0 0.00001 +module_L4b_halfmodule_axial_hole_sensor0 0.00001 +module_L4b_halfmodule_axial_slot_sensor0 0.00001 +module_L4b_halfmodule_stereo_hole_sensor0 0.00001 +module_L4b_halfmodule_stereo_slot_sensor0 0.00001 +module_L5b_halfmodule_axial_hole_sensor0 0.00001 +module_L5b_halfmodule_axial_slot_sensor0 0.00001 +module_L5b_halfmodule_stereo_hole_sensor0 0.00001 +module_L5b_halfmodule_stereo_slot_sensor0 0.00001 +module_L6b_halfmodule_axial_hole_sensor0 0.00001 +module_L6b_halfmodule_axial_slot_sensor0 0.00001 +module_L6b_halfmodule_stereo_hole_sensor0 0.00001 +module_L6b_halfmodule_stereo_slot_sensor0 0.00001 +module_L1t_halfmodule_axial_sensor0 0.00001 +module_L1t_halfmodule_stereo_sensor0 0.00001 +module_L2t_halfmodule_axial_sensor0 0.00001 +module_L2t_halfmodule_stereo_sensor0 0.00001 +module_L3t_halfmodule_axial_sensor0 0.00001 +module_L3t_halfmodule_stereo_sensor0 0.00001 +module_L4t_halfmodule_axial_hole_sensor0 0.00001 +module_L4t_halfmodule_axial_slot_sensor0 0.00001 +module_L4t_halfmodule_stereo_hole_sensor0 0.00001 +module_L4t_halfmodule_stereo_slot_sensor0 0.00001 +module_L5t_halfmodule_axial_hole_sensor0 0.00001 +module_L5t_halfmodule_axial_slot_sensor0 0.00001 +module_L5t_halfmodule_stereo_hole_sensor0 0.00001 +module_L5t_halfmodule_stereo_slot_sensor0 0.00001 +module_L6t_halfmodule_axial_hole_sensor0 0.00001 +module_L6t_halfmodule_axial_slot_sensor0 0.00001 +module_L6t_halfmodule_stereo_hole_sensor0 0.00001 +module_L6t_halfmodule_stereo_slot_sensor0 0.00001 From 61ed3e763b92fa0b916698b71f30d249d63b4d50 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 8 Jan 2025 15:20:54 -0800 Subject: [PATCH 3/4] get rid of all 3d plots --- .../recon/tracking/kalman/KFOutputDriver.java | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java index 083c7186c..c1a7b486a 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java @@ -380,9 +380,9 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda").fill(tanL,eop); aidaKF.histogram2D(eopFolder+"EoP_vs_phi").fill(phi,eop); - aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, - phi, - eop); + // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi").fill(tanL, + // phi, + // eop); if (TriggerModule.inFiducialRegion(cluster)) { @@ -400,9 +400,9 @@ private void doEoPPlots(Track track, Cluster cluster) { aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid").fill(tanL,eop); aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid").fill(phi,eop); - aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(tanL, - phi, - eop); + // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid").fill(tanL, + //phi, + //eop); @@ -534,10 +534,12 @@ private void FillKFTrackPlot(String str, String isTop, String charge, double val aidaKF.histogram2D(str+isTop+charge).fill(valX,valY); } + /* private void FillKFTrackPlot(String str, String isTop, String charge, double valX, double valY, double valZ) { aidaKF.histogram3D(str+isTop).fill(valX,valY,valZ); aidaKF.histogram3D(str+isTop+charge).fill(valX,valY,valZ); } + */ @@ -598,7 +600,7 @@ private void doBasicKFtrack(Track trk, Map sensorHits) //Momentum maps FillKFTrackPlot(trkpFolder+"p_vs_phi",isTop,charge,trackState.getPhi(),trackp); FillKFTrackPlot(trkpFolder+"p_vs_tanLambda",isTop,charge,trackState.getTanLambda(),trackp); - FillKFTrackPlot(trkpFolder+"p_vs_phi_tanLambda",isTop,charge,trackState.getPhi(),trackState.getTanLambda(),trackp); + // FillKFTrackPlot(trkpFolder+"p_vs_phi_tanLambda",isTop,charge,trackState.getPhi(),trackState.getTanLambda(),trackp); double tanLambda = trackState.getTanLambda(); double cosLambda = 1. / (Math.sqrt(1+tanLambda*tanLambda)); @@ -1034,11 +1036,11 @@ private void setupEoPPlots() { aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda",200,-0.1,0.1,200,0,2); aidaKF.histogram2D(eopFolder+"EoP_vs_phi",200,-0.2,0.2,200,0,2); - aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi",200,-0.08,0.08,200,-0.2,0.2,200,0,2); aidaKF.histogram2D(eopFolder+"EoP_vs_tanLambda_fid",200,-0.1,0.1,200,0,2); aidaKF.histogram2D(eopFolder+"EoP_vs_phi_fid",200,-0.2,0.2,200,0,2); - aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid",200,-0.08,0.08,200,-0.2,0.2,200,0,2); + // aidaKF.histogram3D(eopFolder+"EoP_vs_tanLambda_phi_fid",200,-0.08,0.08,200,-0.2,0.2,200,0,2); } @@ -1226,7 +1228,7 @@ private void setupPlots() { aidaKF.histogram2D(trkpFolder+"p_Missing1Hit"+vol+charge,8,0,8,nbins_p,0.0,pmax); aidaKF.histogram2D(trkpFolder+"p_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); aidaKF.histogram2D(trkpFolder+"p_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); - aidaKF.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); + // aidaKF.histogram3D(trkpFolder+"p_vs_phi_tanLambda"+vol+charge, 50,-0.3,0.3,50,-0.2,0.2,100,0.,pmax); aidaKF.histogram2D(trkpFolder+"pT_vs_phi"+vol+charge, nbins_t,-0.3,0.3, nbins_p,0.,pmax); aidaKF.histogram2D(trkpFolder+"pT_vs_tanLambda"+vol+charge,nbins_t,-0.2,0.2,nbins_p,0.,pmax); From 4e243b60950229c2ed4fd9c3fe3552c117944b9a Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Fri, 10 Jan 2025 17:25:35 +0000 Subject: [PATCH 4/4] made comment more useful --- .../recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim index c1d734ef0..0928aa341 100644 --- a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_Killer_Smearer.lcsim @@ -1,10 +1,12 @@