6 modified files
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/SVTHitMCEfficiency.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/SVTHitMCEfficiency.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -2,13 +2,14 @@
import hep.aida.IHistogram1D;
import hep.aida.IHistogram2D;
-import hep.aida.IHistogramFactory;
import hep.aida.IProfile1D;
+import hep.aida.IProfile2D;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.hps.recon.tracking.ShapeFitParameters;
import org.lcsim.detector.tracker.silicon.SiSensor;
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCRelation;
@@ -70,10 +71,16 @@
createLayerPlot(plotDir + "clusterEfficiency", kk, 50, -40, 40.);
createLayerPlot(plotDir + "readoutEfficiency", kk, 50, -40, 40.);
createLayerPlot(plotDir + "rthToClusterEfficiency", kk, 50, -40, 40.);
- createLayerPlot2D(plotDir + "toogoodFits", kk, 200, -100, 100, 100, 0, 20000);
- createLayerPlot2D(plotDir + "goodFits", kk, 200, -100, 100, 100, 0, 20000);
- createLayerPlot2D(plotDir + "badFits", kk, 200, -100, 100, 100, 0, 20000);
+ createLayerPlot2D(plotDir + "clusterEfficiency2D", kk, 50, -40, 40., 16, 0.5, 16.5);
+ createLayerPlot2D(plotDir + "rthToClusterEfficiency2D", kk, 50, -40, 40., 16, 0.5, 16.5);
+ createLayerPlot2D(plotDir + "allFits", kk, 200, -100, 100, 100, 0, 20000);
+// createLayerPlot2D(plotDir + "toogoodFits", kk, 200, -100, 100, 100, 0, 20000);
+// createLayerPlot2D(plotDir + "goodFits", kk, 200, -100, 100, 100, 0, 20000);
+// createLayerPlot2D(plotDir + "badFits", kk, 200, -100, 100, 100, 0, 20000);
+ createLayerPlot2D(plotDir + "fitT0ChiProb", kk, 200, -100, 100, 100, 0, 1.0);
+ createLayerPlot2D(plotDir + "fitAmpChiProb", kk, 200, 0, 20000, 100, 0, 1.0);
createLayerPlot1D(plotDir + "signalClusterT0", kk, 500, -100, 100);
+ createLayerPlot2D(plotDir + "badClusterFits", kk, 200, -100, 100, 100, 0, 20000);
}
resetEfficiencyMap();
}
@@ -141,9 +148,14 @@
}
//relational tables from raw and fitted tracker hits to sim hit
+ RelationalTable rthtofit = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED);
RelationalTable fittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<FittedRawTrackerHit> fittedTrackerHits = event.get(FittedRawTrackerHit.class, fittedTrackerHitCollectionName);
for (FittedRawTrackerHit hit : fittedTrackerHits) {
+ ShapeFitParameters oldfit = (ShapeFitParameters) rthtofit.to(hit.getRawTrackerHit());
+ if (oldfit == null || Math.abs(oldfit.getT0()) > Math.abs(hit.getT0())) {
+ rthtofit.add(hit.getRawTrackerHit(), hit.getShapeFitParameters());
+ }
Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(hit.getRawTrackerHit());
for (SimTrackerHit simhit : simTrackerHits) {
fittomc.add(hit, simhit);
@@ -160,24 +172,39 @@
}
if (signalHit != null) {
// System.out.format("chiprob %f, t0 %f, A %f\n", signalHit.getShapeFitParameters().getChiProb(), signalHit.getT0(), signalHit.getAmp());
- if (signalHit.getShapeFitParameters().getChiProb() > 0.95) {
- getLayerPlot2D(plotDir + "toogoodFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
- } else if (signalHit.getShapeFitParameters().getChiProb() < 0.05) {
- getLayerPlot2D(plotDir + "badFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
- } else {
- getLayerPlot2D(plotDir + "goodFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
- }
+ getLayerPlot2D(plotDir + "allFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
+ getLayerPlot2D(plotDir + "fitT0ChiProb", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getShapeFitParameters().getChiProb());
+ getLayerPlot2D(plotDir + "fitAmpChiProb", simhit.getLayer()).fill(signalHit.getAmp(), signalHit.getShapeFitParameters().getChiProb());
+// if (signalHit.getShapeFitParameters().getChiProb() > 0.95) {
+// getLayerPlot2D(plotDir + "toogoodFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
+// } else if (signalHit.getShapeFitParameters().getChiProb() < 0.05) {
+// getLayerPlot2D(plotDir + "badFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
+// } else {
+// getLayerPlot2D(plotDir + "goodFits", simhit.getLayer()).fill(signalHit.getT0(), signalHit.getAmp());
+// }
}
int gotCluster = 0;
+ int[] gotClusterAtTime = new int[16];
Set<SiTrackerHitStrip1D> clusters = clustertosimhit.allTo(simhit);
if (clusters != null) {
for (SiTrackerHitStrip1D clust : clusters) {
getLayerPlot1D(plotDir + "signalClusterT0", simhit.getLayer()).fill(clust.getTime());
+ for (int i = 0; i < 16; i++) {
+ if (Math.abs(clust.getTime()) < i + 1) {
+ gotClusterAtTime[i] = 1;
+ }
+ }
if (Math.abs(clust.getTime()) < t0Cut) {
gotCluster = 1;
+ } else {
+ for (RawTrackerHit rth : clust.getRawHits()) {
+ ShapeFitParameters fit = (ShapeFitParameters) rthtofit.to(rth);
+ getLayerPlot2D(plotDir + "badClusterFits", simhit.getLayer()).fill(fit.getT0(), fit.getAmp());
+ }
}
+
}
}
Set<RawTrackerHit> rawhits = rawtomc.allTo(simhit);
@@ -189,11 +216,21 @@
if (gotRawHit == 1) {
getLayerPlot(plotDir + "rthToClusterEfficiency", simhit.getLayer()).fill(y, gotCluster);
}
+ for (int i = 0; i < 16; i++) {
+ getLayerPlot2D(plotDir + "clusterEfficiency2D", simhit.getLayer()).fill(y, i + 1, gotClusterAtTime[i]);
+ if (gotRawHit == 1) {
+ getLayerPlot2D(plotDir + "rthToClusterEfficiency2D", simhit.getLayer()).fill(y, i + 1, gotClusterAtTime[i]);
+ }
+ }
}
}
@Override
public void fillEndOfRunPlots() {
+ for (int kk = 1; kk < 13; kk++) {
+ getMean2D(getLayerPlot2D(plotDir + "clusterEfficiency2D", kk));
+ getMean2D(getLayerPlot2D(plotDir + "rthToClusterEfficiency2D", kk));
+ }
}
@Override
@@ -224,6 +261,33 @@
return aida.profile1D(prefix + "_layer" + layer, nchan, min, max);
}
+ private void getMean2D(IHistogram2D hist2D) {
+ int nx = hist2D.xAxis().bins();
+ int ny = hist2D.yAxis().bins();
+ double[][] means = new double[nx][ny];
+ for (int ix = 0; ix < nx; ix++) {
+ for (int iy = 0; iy < ny; iy++) {
+ means[ix][iy] = hist2D.binHeight(ix, iy) / hist2D.binEntries(ix, iy);
+ }
+ }
+ hist2D.reset();
+ for (int ix = 0; ix < nx; ix++) {
+ for (int iy = 0; iy < ny; iy++) {
+ double x = hist2D.xAxis().binCenter(ix);
+ double y = hist2D.yAxis().binCenter(iy);
+ hist2D.fill(x, y, means[ix][iy]);
+ }
+ }
+ }
+
+ private IProfile2D getLayerProfile2D(String prefix, int layer) {
+ return aida.profile2D(prefix + "_layer" + layer);
+ }
+
+ private IProfile2D createLayerProfile2D(String prefix, int layer, int nx, double minX, double maxX, int ny, double minY, double maxY) {
+ return aida.profile2D(prefix + "_layer" + layer, nx, minX, maxX, ny, minY, maxY);
+ }
+
private IHistogram1D getLayerPlot1D(String prefix, int layer) {
return aida.histogram1D(prefix + "_layer" + layer);
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/DataTrackerHitDriver.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/DataTrackerHitDriver.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -38,6 +38,9 @@
private double clusterSeedThreshold = 4.0;
private double clusterNeighborThreshold = 3.0;
private double clusterThreshold = 4.0;
+ private double meanTime = 24.0;
+ private double timeWindow = 48.0;
+ private double neighborDeltaT = 24.0;
private int clusterMaxSize = 10;
private int clusterCentralStripAveragingThreshold = 4;
// Clustering errors by number of TrackerHits.
@@ -85,6 +88,18 @@
this.clusterThreshold = clusterThreshold;
}
+ public void setMeanTime(double meanTime) {
+ this.meanTime = meanTime;
+ }
+
+ public void setTimeWindow(double timeWindow) {
+ this.timeWindow = timeWindow;
+ }
+
+ public void setNeighborDeltaT(double neighborDeltaT) {
+ this.neighborDeltaT = neighborDeltaT;
+ }
+
public void setClusterMaxSize(int clusterMaxSize) {
this.clusterMaxSize = clusterMaxSize;
}
@@ -154,6 +169,9 @@
stripClusteringAlgo.setSeedThreshold(clusterSeedThreshold);
stripClusteringAlgo.setNeighborThreshold(clusterNeighborThreshold);
stripClusteringAlgo.setClusterThreshold(clusterThreshold);
+ stripClusteringAlgo.setMeanTime(meanTime);
+ stripClusteringAlgo.setTimeWindow(timeWindow);
+ stripClusteringAlgo.setNeighborDeltaT(neighborDeltaT);
// hitMaker=new HPSFittedRawTrackerHitMaker(shaperFit);
// Create the clusterers and set hit-making parameters.
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -1,6 +1,7 @@
package org.hps.recon.tracking;
import java.util.ArrayList;
+import java.util.Collection;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedList;
@@ -28,6 +29,7 @@
private double _cluster_threshold;
private double _meanTime = 24;
private double _timeWindow = 48;
+ private double _neighborDeltaT = Double.POSITIVE_INFINITY;
private final double _minChiProb = Gamma.regularizedGammaQ(4, 20);
/**
@@ -47,6 +49,18 @@
_cluster_threshold = cluster_threshold;
}
+ public void setMeanTime(double _meanTime) {
+ this._meanTime = _meanTime;
+ }
+
+ public void setTimeWindow(double _timeWindow) {
+ this._timeWindow = _timeWindow;
+ }
+
+ public void setNeighborDeltaT(double _neighborDeltaT) {
+ this._neighborDeltaT = _neighborDeltaT;
+ }
+
/**
* Instantiate NearestNeighborRMS with default thresholds:
*
@@ -101,8 +115,10 @@
// Create maps that show the channel status and relate the channel number to the raw hit
// and vice versa
int mapsize = 2 * base_hits.size();
- Map<Integer, Boolean> clusterable = new HashMap<Integer, Boolean>(mapsize);
- Map<FittedRawTrackerHit, Integer> hit_to_channel = new HashMap<FittedRawTrackerHit, Integer>(mapsize);
+// Map<Integer, Boolean> clusterable = new HashMap<Integer, Boolean>(mapsize);
+ Set<Integer> clusterableSet = new HashSet<Integer>(mapsize);
+
+// Map<FittedRawTrackerHit, Integer> hit_to_channel = new HashMap<FittedRawTrackerHit, Integer>(mapsize);
Map<Integer, FittedRawTrackerHit> channel_to_hit = new HashMap<Integer, FittedRawTrackerHit>(mapsize);
// Create list of channel numbers to be used as cluster seeds
@@ -119,11 +135,10 @@
IIdentifier id = rth.getIdentifier();
int channel_number = sid_helper.getElectrodeValue(id);
- // Check for duplicate RawTrackerHit
- if (hit_to_channel.containsKey(base_hit)) {
- throw new RuntimeException("Duplicate hit: " + id.toString());
- }
-
+// // Check for duplicate RawTrackerHit
+// if (hit_to_channel.containsKey(base_hit)) {
+// throw new RuntimeException("Duplicate hit: " + id.toString());
+// }
// Check for duplicate RawTrackerHits or channel numbers
if (channel_to_hit.containsKey(channel_number)) {
// throw new RuntimeException("Duplicate channel number: "+channel_number);
@@ -136,15 +151,17 @@
}
// Add this hit to the maps that relate channels and hits
- hit_to_channel.put(base_hit, channel_number);
+// hit_to_channel.put(base_hit, channel_number);
channel_to_hit.put(channel_number, base_hit);
// Get the signal from the readout chip
double signal = base_hit.getAmp();
double noiseRMS = HPSSVTCalibrationConstants.getNoise((SiSensor) rth.getDetectorElement(), channel_number);
- double time = base_hit.getT0();
+
// Mark this hit as available for clustering if it is above the neighbor threshold
- clusterable.put(channel_number, signal / noiseRMS >= _neighbor_threshold);
+ if (signal / noiseRMS >= _neighbor_threshold && passChisqCut(base_hit)) {
+ clusterableSet.add(channel_number);
+ }
// Add this hit to the list of seeds if it is above the seed threshold
if (signal / noiseRMS >= _seed_threshold && passTimingCut(base_hit) && passChisqCut(base_hit)) {
@@ -159,7 +176,10 @@
for (int seed_channel : cluster_seeds) {
// First check if this hit is still available for clustering
- if (!clusterable.get(seed_channel)) {
+// if (!clusterable.get(seed_channel)) {
+// continue;
+// }
+ if (!clusterableSet.contains(seed_channel)) {
continue;
}
@@ -167,13 +187,15 @@
List<FittedRawTrackerHit> cluster = new ArrayList<FittedRawTrackerHit>();
double cluster_signal = 0.;
double cluster_noise_squared = 0.;
+ double cluster_weighted_time = 0.;
// Create a queue to hold channels whose neighbors need to be checked for inclusion
LinkedList<Integer> unchecked = new LinkedList<Integer>();
// Add the seed channel to the unchecked list and mark it as unavailable for clustering
unchecked.addLast(seed_channel);
- clusterable.put(seed_channel, false);
+// clusterable.put(seed_channel, false);
+ clusterableSet.remove(seed_channel);
// Check the neighbors of channels added to the cluster
while (unchecked.size() > 0) {
@@ -181,34 +203,34 @@
// Pull the next channel off the queue and add it's hit to the cluster
int clustered_cell = unchecked.removeFirst();
cluster.add(channel_to_hit.get(clustered_cell));
- cluster_signal += channel_to_hit.get(clustered_cell).getAmp();
- cluster_noise_squared += Math.pow(HPSSVTCalibrationConstants.getNoise((SiSensor) (channel_to_hit.get(clustered_cell)).getRawTrackerHit().getDetectorElement(), clustered_cell), 2);
+ FittedRawTrackerHit hit = channel_to_hit.get(clustered_cell);
+ cluster_signal += hit.getAmp();
+ cluster_noise_squared += Math.pow(HPSSVTCalibrationConstants.getNoise((SiSensor) hit.getRawTrackerHit().getDetectorElement(), clustered_cell), 2);
+ cluster_weighted_time += hit.getT0() * hit.getAmp();
// cluster_noise_squared +=0; //need to get the noise from the calib. const. class
// Get the neigbor channels
// Set<Integer> neighbor_channels =
// electrodes.getNearestNeighborCells(clustered_cell);
- Set<Integer> neighbor_channels = getNearestNeighborCells(clustered_cell);
+ Collection<Integer> neighbor_channels = getNearestNeighborCells(clustered_cell);
// Now loop over the neighbors and see if we can add them to the cluster
for (int channel : neighbor_channels) {
- // Get the status of this channel
- Boolean addhit = clusterable.get(channel);
-
- // If the map entry is null, there is no raw hit for this channel
- if (addhit == null) {
+ // Check if this neighbor channel is still available for clustering
+ if (!clusterableSet.contains(channel)) {
continue;
}
- // Check if this neighbor channel is still available for clustering
- if (!addhit) {
+ FittedRawTrackerHit neighbor_hit = channel_to_hit.get(channel);
+ if (Math.abs(neighbor_hit.getT0() - cluster_weighted_time / cluster_signal) > _neighborDeltaT) {
+// System.out.format("new hit t0 %f, cluster t0 %f\n", neighbor_hit.getT0(), cluster_weighted_time / cluster_signal);
continue;
}
// Add channel to the list of unchecked clustered channels
// and mark it unavailable for clustering
unchecked.addLast(channel);
- clusterable.put(channel, false);
+ clusterableSet.remove(channel);
} // end of loop over neighbor cells
} // end of loop over unchecked cells
@@ -226,33 +248,26 @@
}
private boolean passTimingCut(FittedRawTrackerHit hit) {
-
- boolean pass = false;
double time = hit.getT0();
- if (Math.abs(time - _meanTime) < _timeWindow) {
- pass = true;
- }
-
- return pass;
+ return (Math.abs(time - _meanTime) < _timeWindow);
}
private boolean passChisqCut(FittedRawTrackerHit hit) {
return hit.getShapeFitParameters().getChiProb() > _minChiProb;
}
- public int getNeighborCell(int cell, int ncells_0, int ncells_1) {
- int neighbor_cell = cell + ncells_0;
- if (isValidCell(neighbor_cell)) {
- return neighbor_cell;
- } else {
- return -1;
- }
- }
-
- public Set<Integer> getNearestNeighborCells(int cell) {
- Set<Integer> neighbors = new HashSet<Integer>();
+// public int getNeighborCell(int cell, int ncells_0, int ncells_1) {
+// int neighbor_cell = cell + ncells_0;
+// if (isValidCell(neighbor_cell)) {
+// return neighbor_cell;
+// } else {
+// return -1;
+// }
+// }
+ public Collection<Integer> getNearestNeighborCells(int cell) {
+ Collection<Integer> neighbors = new ArrayList<Integer>(2);
for (int ineigh = -1; ineigh <= 1; ineigh = ineigh + 2) {
- int neighbor_cell = getNeighborCell(cell, ineigh, 0);
+ int neighbor_cell = cell + ineigh;
if (isValidCell(neighbor_cell)) {
neighbors.add(neighbor_cell);
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -10,6 +10,7 @@
import org.lcsim.detector.tracker.silicon.SiSensor;
import org.lcsim.event.EventHeader;
import org.lcsim.event.RawTrackerHit;
+import org.lcsim.geometry.Detector;
import org.lcsim.lcio.LCIOConstants;
import org.lcsim.recon.cat.util.Const;
import org.lcsim.util.Driver;
@@ -85,7 +86,7 @@
}
@Override
- public void startOfData() {
+ public void detectorChanged(Detector detector) {
_shaper.setDebug(debug);
if (rawHitCollectionName == null) {
throw new RuntimeException("The parameter ecalCollectionName was not set!");
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -46,10 +46,10 @@
this.nPulses = nPulses;
amplitudes = new double[nPulses];
amplitudeErrors = new double[nPulses];
- System.setErr(new PrintStream(new OutputStream() {
- public void write(int b) {
- }
- }));
+// System.setErr(new PrintStream(new OutputStream() {
+// public void write(int b) {
+// }
+// }));
}
public void setDebug(boolean debug) {
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java 2014-09-03 18:35:04 UTC (rev 946)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java 2014-09-03 21:00:22 UTC (rev 947)
@@ -15,6 +15,9 @@
ShaperLinearFitAlgorithm twoPulseFitter = new ShaperLinearFitAlgorithm(2);
private boolean debug = false;
private double refitThreshold = 0.5;
+ private int totalFits = 0;
+ private int refitAttempts = 0;
+ private int refitsAccepted = 0;
public ShaperPileupFitAlgorithm() {
}
@@ -26,18 +29,26 @@
public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, HPSSVTCalibrationConstants.ChannelConstants constants) {
Collection<ShapeFitParameters> fittedPulses = onePulseFitter.fitShape(rth, constants);
double singlePulseChiProb = fittedPulses.iterator().next().getChiProb();
+ totalFits++;
if (singlePulseChiProb < refitThreshold) {
+ refitAttempts++;
Collection<ShapeFitParameters> doublePulse = twoPulseFitter.fitShape(rth, constants);
double doublePulseChiProb = doublePulse.iterator().next().getChiProb();
if (doublePulseChiProb > singlePulseChiProb) {
+ refitsAccepted++;
fittedPulses = doublePulse;
}
}
+ if (debug && totalFits % 10000 == 0) {
+ System.out.format("%d fits, %d refit attempts, %d refits accepted\n", totalFits, refitAttempts, refitsAccepted);
+ }
return fittedPulses;
}
public void setDebug(boolean debug) {
this.debug = debug;
+ onePulseFitter.setDebug(debug);
+ twoPulseFitter.setDebug(debug);
}
}
SVNspam 0.1