Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking on MAIN | |||
HPSShaperAnalyticFitAlgorithm.java | +86 | added 1.1 | |
HPSShapeFitParameters.java | +15 | -3 | 1.2 -> 1.3 |
HPSSVTCalibrationConstants.java | +131 | -88 | 1.1 -> 1.2 |
HPSRawTrackerHitFitterDriver.java | +6 | -4 | 1.2 -> 1.3 |
HPSSVTConstants.java | +3 | -1 | 1.2 -> 1.3 |
HPSShaperFitAlgorithm.java | +2 | -1 | 1.2 -> 1.3 |
HPSTrackerHitMaker.java | -41 | 1.2 removed | |
+243 | -138 |
add one fit algorithm
diff -N HPSShaperAnalyticFitAlgorithm.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HPSShaperAnalyticFitAlgorithm.java 24 Apr 2012 23:32:31 -0000 1.1 @@ -0,0 +1,86 @@
+package org.lcsim.hps.recon.tracking; + +import org.lcsim.event.RawTrackerHit; +import org.lcsim.hps.recon.tracking.HPSSVTCalibrationConstants.ChannelConstants; + +/** + * Fast fitter; currently only fits single hits. Uses Tp from ChannelConstants; fits values and errors for T0 and amplitude. + * @author meeg + * @version $Id: HPSShaperAnalyticFitAlgorithm.java,v 1.1 2012/04/24 23:32:31 meeg Exp $ + */ +public class HPSShaperAnalyticFitAlgorithm implements HPSShaperFitAlgorithm { + + public HPSShapeFitParameters fitShape(RawTrackerHit rth, ChannelConstants constants) { + double minChisq = Double.POSITIVE_INFINITY; + int bestStart = 0; + HPSShapeFitParameters fit = new HPSShapeFitParameters(); + for (int i = 0; i < rth.getADCValues().length - 2; i++) { + double chisq = fitSection(rth, constants, fit, i); + if (chisq < minChisq) { + minChisq = chisq / (rth.getADCValues().length - 2 - i); //TODO: p-value would be better here + bestStart = i; + } + } + fitSection(rth, constants, fit, bestStart); + return fit; + } + + private double fitSection(RawTrackerHit rth, ChannelConstants constants, HPSShapeFitParameters fit, int start) { + int length = rth.getADCValues().length - start + 1; + double[] y = new double[length]; + double[] t = new double[length]; + + for (int i = 0; i < length; i++) { + y[i] = rth.getADCValues()[start + i] - constants.getPedestal(); + t[i] = HPSSVTConstants.SAMPLE_INTERVAL * i; + } + + double[] p = new double[length]; + double[] a = new double[length]; + for (int i = 0; i < length; i++) { + p[i] = y[i] / constants.getNoise(); + a[i] = Math.exp(1 - t[i] / constants.getTp()) / (constants.getTp() * constants.getNoise()); + } + + double pa, aatt, pat, aat, aa; + pa = 0; + aatt = 0; + pat = 0; + aat = 0; + aa = 0; + for (int i = 0; i < length; i++) { + pa += p[i] * a[i]; + aatt += a[i] * a[i] * t[i] * t[i]; + pat += p[i] * t[i] * t[i]; + aat += a[i] * a[i] * t[i]; + aa += a[i] * a[i]; + } + + double t0 = (pa * aatt - pat * aat) / (pa * aat - aa * pat); + double A = pa / (Math.exp(t0 / constants.getTp()) * (aat - t0 * aa)); + + double time_var = 0; + double height_var = 0; + for (int i = 0; i < length; i++) { + double dt_dp = a[i] * (aatt - t[i] * aat - t0 * (aat - t[i] * aa)) / (pa * aat - aa * pat); + double dh_dp = (a[i] * Math.exp(-1.0 * t0 / constants.getTp()) + A * dt_dp * aa) / (aat - t0 * aa) - A * dt_dp / constants.getTp(); + time_var += dt_dp * dt_dp; + height_var += dh_dp * dh_dp; +// covar += dt_dp*dh_dp; + } + + fit.setAmp(A); + fit.setAmpErr(Math.sqrt(height_var)); + fit.setT0(t0); + fit.setT0Err(Math.sqrt(time_var)); + fit.setTp(constants.getTp()); + + double chisq = 0; + for (int i = 0; i < length; i++) { + double fit_y = A * (t[i] / constants.getTp()) * Math.exp(1 - t[i] / constants.getTp()); + chisq += Math.pow((fit_y - y[i]) / constants.getNoise(), 2); + } + + return chisq; + } +}
diff -u -r1.2 -r1.3 --- HPSShapeFitParameters.java 24 Apr 2012 16:23:38 -0000 1.2 +++ HPSShapeFitParameters.java 24 Apr 2012 23:32:31 -0000 1.3 @@ -35,10 +35,22 @@
_amp = amp; }
- public void setTP(double tp) {
+ public void setTp(double tp) {
_tp = tp; }
+ public void setAmpErr(double _ampErr) { + this._ampErr = _ampErr; + } + + public void setT0Err(double _t0Err) { + this._t0Err = _t0Err; + } + + public void setTpErr(double _tpErr) { + this._tpErr = _tpErr; + } +
public double getT0() { return _t0; }
@@ -47,7 +59,7 @@
return _amp; }
- public double getTP() {
+ public double getTp() {
return _tp; }
@@ -59,7 +71,7 @@
return _ampErr; }
- public double getTPErr() {
+ public double getTpErr() {
return _tpErr; }
diff -u -r1.1 -r1.2 --- HPSSVTCalibrationConstants.java 24 Apr 2012 14:06:41 -0000 1.1 +++ HPSSVTCalibrationConstants.java 24 Apr 2012 23:32:31 -0000 1.2 @@ -16,96 +16,139 @@
* * @author mgraham */
-public class HPSSVTCalibrationConstants{
+public class HPSSVTCalibrationConstants {
- // Map from Sensor to Hybrid/FPGA pair - public static Map<Pair<SiSensor /* - * Sensor - */, Integer> /* - * Channel Number - */, Double /* - * constant - */> noise; - public static Map<Pair<SiSensor /* - * Sensor - */, Integer> /* - * Channel Number - */, Double /* - * constant - */> pedestal; - public static Map<Pair<SiSensor /* - * Sensor - */, Integer> /* - * Channel Number - */, Double /* - * constant - */> tShaping; - - /** - * Default Constructor. - */ - public HPSSVTCalibrationConstants() { - pedestal = new HashMap<Pair<SiSensor, Integer>, Double>(); - noise = new HashMap<Pair<SiSensor, Integer>, Double>(); - tShaping = new HashMap<Pair<SiSensor, Integer>, Double>();
+ // Map from Sensor to Hybrid/FPGA pair + public static Map<Pair<SiSensor /* + * Sensor + */, Integer> /* + * Channel Number + */, Double /* + * constant + */> noise; + public static Map<Pair<SiSensor /* + * Sensor + */, Integer> /* + * Channel Number + */, Double /* + * constant + */> pedestal; + public static Map<Pair<SiSensor /* + * Sensor + */, Integer> /* + * Channel Number + */, Double /* + * constant + */> tShaping; + + /** + * Default Constructor. + */ + public HPSSVTCalibrationConstants() { + pedestal = new HashMap<Pair<SiSensor, Integer>, Double>(); + noise = new HashMap<Pair<SiSensor, Integer>, Double>(); + tShaping = new HashMap<Pair<SiSensor, Integer>, Double>();
- }
+ }
- public void loadCalibrationConstants(String calibFile) { - //write something here to read in constants from calibration fiel
+ public void loadCalibrationConstants(String calibFile) { + //write something here to read in constants from calibration fiel
/*
- ConditionsManager conditions = ConditionsManager.defaultInstance(); - StreamTokenizer tok = null; - try { - tok = new StreamTokenizer(conditions.getRawConditions(calibFile).getReader()); - } catch (IOException e) { - throw new RuntimeException("couldn't get DAQ map from conditions manager", e); - } - tok.commentChar('#'); - tok.eolIsSignificant(false); - tok.parseNumbers(); - try { - while (tok.nextToken() != StreamTokenizer.TT_EOF) { - tok.pushBack(); - int x = (int) getNumber(tok); - int y = (int) getNumber(tok); - int crate = (int) getNumber(tok); - short slot = (short) getNumber(tok); - short channel = (short) getNumber(tok); - addMapEntry(enc.getID(), getDaqID(crate, slot, channel)); - pedestal.put(sensorChannel, physicalID); - noise.put(sensorChannel, daqID); - // System.out.printf("x = %d, y = %d, crate = %d,slot = % d, channel = % d\n", x, y, crate, slot, channel); // - //System.out.printf("physicalID = %d, daqID = %d\n", enc.getID(), getDaqID(crate, slot, channel)); - } - } catch (IOException e) { - throw new RuntimeException("couldn't parse DAQ map", e); - } - * - */ - - } - - public static double getNoise(SiSensor sensor, int channel) { - Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); - return noise.get(sensorChannel); - } - - public static double getPedestal(SiSensor sensor, int channel) { - Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); - return pedestal.get(sensorChannel); - } - - public static double getTShaping(SiSensor sensor, int channel) { - Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); - return tShaping.get(sensorChannel); - } - - private double getNumber(StreamTokenizer tok) throws IOException { - if (tok.nextToken() == StreamTokenizer.TT_NUMBER) { - return tok.nval; - } else { - throw new IOException("expected an int from DAQ map, instead got " + tok); - } - }
+ * ConditionsManager conditions = + * ConditionsManager.defaultInstance(); StreamTokenizer tok = + * null; try { tok = new + * StreamTokenizer(conditions.getRawConditions(calibFile).getReader()); + * } catch (IOException e) { throw new + * RuntimeException("couldn't get DAQ map from conditions + * manager", e); } tok.commentChar('#'); + * tok.eolIsSignificant(false); tok.parseNumbers(); try { while + * (tok.nextToken() != StreamTokenizer.TT_EOF) { tok.pushBack(); + * int x = (int) getNumber(tok); int y = (int) getNumber(tok); + * int crate = (int) getNumber(tok); short slot = (short) + * getNumber(tok); short channel = (short) getNumber(tok); + * addMapEntry(enc.getID(), getDaqID(crate, slot, channel)); + * pedestal.put(sensorChannel, physicalID); + * noise.put(sensorChannel, daqID); // System.out.printf("x = + * %d, y = %d, crate = %d,slot = % d, channel = % d\n", x, y, + * crate, slot, channel); // //System.out.printf("physicalID = + * %d, daqID = %d\n", enc.getID(), getDaqID(crate, slot, + * channel)); } } catch (IOException e) { throw new + * RuntimeException("couldn't parse DAQ map", e); } + * + */ + } + + public static double getNoise(SiSensor sensor, int channel) { + Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); + return noise.get(sensorChannel); + } + + public static double getPedestal(SiSensor sensor, int channel) { + Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); + return pedestal.get(sensorChannel); + } + + public static double getTShaping(SiSensor sensor, int channel) { + Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); + return tShaping.get(sensorChannel); + } + + private double getNumber(StreamTokenizer tok) throws IOException { + if (tok.nextToken() == StreamTokenizer.TT_NUMBER) { + return tok.nval; + } else { + throw new IOException("expected an int from DAQ map, instead got " + tok); + } + } + + public static ChannelConstants getChannelConstants(SiSensor sensor, int channel) { + Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel); + ChannelConstants constants = new ChannelConstants(); + constants.setNoise(noise.get(sensorChannel)); + constants.setPedestal(pedestal.get(sensorChannel)); + constants.setTp(tShaping.get(sensorChannel)); + + return constants; + } + + //class to hold calibration constants for a channel; leave fields NaN or null if not known + public static class ChannelConstants { + + private double pedestal = Double.NaN; + private double tp = Double.NaN; + private double noise = Double.NaN; + private double[][] pulseShape = null; + + public double getNoise() { + return noise; + } + + public void setNoise(double noise) { + this.noise = noise; + } + + public double getPedestal() { + return pedestal; + } + + public void setPedestal(double pedestal) { + this.pedestal = pedestal; + } + + public double[][] getPulseShape() { + return pulseShape; + } + + public void setPulseShape(double[][] pulseShape) { + this.pulseShape = pulseShape; + } + + public double getTp() { + return tp; + } + + public void setTp(double tp) { + this.tp = tp; + } + }
}
diff -u -r1.2 -r1.3 --- HPSRawTrackerHitFitterDriver.java 24 Apr 2012 18:01:49 -0000 1.2 +++ HPSRawTrackerHitFitterDriver.java 24 Apr 2012 23:32:31 -0000 1.3 @@ -6,8 +6,10 @@
import java.util.ArrayList; import java.util.List;
+import org.lcsim.detector.tracker.silicon.SiSensor;
import org.lcsim.event.EventHeader; import org.lcsim.event.RawTrackerHit;
+import org.lcsim.hps.recon.tracking.HPSSVTCalibrationConstants.ChannelConstants;
import org.lcsim.util.Driver; import org.lcsim.util.lcio.LCIOConstants;
@@ -25,9 +27,8 @@
private int relationFlags = 0; public void setFitAlgorithm(String fitAlgorithm) {
- //TODO: uncomment once this actually exists -// if (fitAlgorithm=="Analytic") -// _shaper = new HPSShaperAnalyticFitter();
+ if (fitAlgorithm=="Analytic") + _shaper = new HPSShaperAnalyticFitAlgorithm();
} public void setFitCollectionName(String fitCollectionName) {
@@ -60,7 +61,8 @@
// Make a fitted hit from this cluster for (RawTrackerHit hit : rawHits) {
- HPSShapeFitParameters fit = _shaper.fitShape(hit);
+ ChannelConstants constants = HPSSVTCalibrationConstants.getChannelConstants((SiSensor) hit.getDetectorElement(), hit.getIdentifierFieldValue("strip")); + HPSShapeFitParameters fit = _shaper.fitShape(hit,constants);
fits.add(fit); HPSFittedRawTrackerHit hth = new HPSFittedRawTrackerHit(hit, fit); hits.add(hth);
diff -u -r1.2 -r1.3 --- HPSSVTConstants.java 22 Apr 2012 21:04:11 -0000 1.2 +++ HPSSVTConstants.java 24 Apr 2012 23:32:31 -0000 1.3 @@ -3,7 +3,7 @@
/** * * @author Omar Moreno <[log in to unmask]>
- * @version $Id: HPSSVTConstants.java,v 1.2 2012/04/22 21:04:11 omoreno Exp $
+ * @version $Id: HPSSVTConstants.java,v 1.3 2012/04/24 23:32:31 meeg Exp $
* */
@@ -20,6 +20,8 @@
// Total number of shaper signal samples obtained public static final int TOTAL_NUMBER_OF_SAMPLES = 6;
+ public static final double SAMPLE_INTERVAL = 24.0; +
// MASK used to encode and decode the SVT data public static final int TEMP_MASK = 0xFFFF; public static final int HYBRID_MASK = 0x3;
diff -u -r1.2 -r1.3 --- HPSShaperFitAlgorithm.java 24 Apr 2012 16:06:25 -0000 1.2 +++ HPSShaperFitAlgorithm.java 24 Apr 2012 23:32:31 -0000 1.3 @@ -5,6 +5,7 @@
package org.lcsim.hps.recon.tracking; import org.lcsim.event.RawTrackerHit;
+import org.lcsim.hps.recon.tracking.HPSSVTCalibrationConstants.ChannelConstants;
/** *
@@ -12,6 +13,6 @@
*/ public interface HPSShaperFitAlgorithm {
- public HPSShapeFitParameters fitShape(RawTrackerHit rth);
+ public HPSShapeFitParameters fitShape(RawTrackerHit rth, ChannelConstants constants);
}
diff -N HPSTrackerHitMaker.java --- HPSTrackerHitMaker.java 24 Apr 2012 16:23:38 -0000 1.2 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,41 +0,0 @@
-/* - * To change this template, choose Tools | Templates - * and open the template in the editor. - */ -package org.lcsim.hps.recon.tracking; - -import java.util.ArrayList; -import java.util.List; -import org.lcsim.event.RawTrackerHit; - -/** - * - * @author mgraham - */ -public class HPSTrackerHitMaker { - - private static String _NAME = "HPSTrackerHitMaker"; - private HPSShaperFitAlgorithm _shaper; - - /** - * Creates a new instance of RawTrackerHitMaker - */ - public HPSTrackerHitMaker(HPSShaperFitAlgorithm shaperFit) { - _shaper = shaperFit; - } - - public List<HPSFittedRawTrackerHit> makeHits(List<RawTrackerHit> raw_hits) { - List<HPSFittedRawTrackerHit> hits = new ArrayList<HPSFittedRawTrackerHit>(); - // Make a pixel hit from this cluster - for (RawTrackerHit hit : raw_hits) { - HPSShapeFitParameters fit = _shaper.fitShape(hit); - HPSFittedRawTrackerHit hth = new HPSFittedRawTrackerHit(hit, fit); - hits.add(hth); - } - return hits; - } - - public String getName() { - return _NAME; - } -}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1