Commit in hps-java/src/main on MAIN | |||
java/org/lcsim/hps/recon/tracking/HPSShaperAnalyticFitAlgorithm.java | +83 | -73 | 1.4 -> 1.5 |
/HPSShapeFitParameters.java | +129 | -123 | 1.4 -> 1.5 |
/HPSSVTCalibrationConstants.java | +16 | -6 | 1.4 -> 1.5 |
/TrackerReconDriver.java | +2 | -1 | 1.15 -> 1.16 |
/HPSRawTrackerHitFitterDriver.java | +56 | -54 | 1.6 -> 1.7 |
/HPSHelicalTrackHitDriver.java | +4 | -11 | 1.4 -> 1.5 |
resources/org/lcsim/hps/steering/OnlineTracking.lcsim | +4 | -1 | 1.2 -> 1.3 |
+294 | -269 |
fix AnalyticFitAlgorithm; works now
diff -u -r1.4 -r1.5 --- HPSShaperAnalyticFitAlgorithm.java 25 Apr 2012 18:01:32 -0000 1.4 +++ HPSShaperAnalyticFitAlgorithm.java 28 Apr 2012 19:35:09 -0000 1.5 @@ -4,83 +4,93 @@
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.
+ * 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.4 2012/04/25 18:01:32 mgraham Exp $
+ * @version $Id: HPSShaperAnalyticFitAlgorithm.java,v 1.4 2012/04/25 18:01:32 + * mgraham 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; - bestStart = i; - } - } - fit.setChiSq(fitSection(rth, constants, fit, bestStart)); - return fit; - } - - private double fitSection(RawTrackerHit rth, ChannelConstants constants, HPSShapeFitParameters fit, int start) { - int length = rth.getADCValues().length - start; - 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;
+ 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; + bestStart = i; + } + } + fit.setChiSq(fitSection(rth, constants, fit, bestStart)); + System.out.println("start=" + bestStart + ", fit " + fit); + return fit; + } + + private double fitSection(RawTrackerHit rth, ChannelConstants constants, HPSShapeFitParameters fit, int start) { + int length = rth.getADCValues().length - start; + 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] * a[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 / (length - 2); //TODO: p-value would be better here - }
+ 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 < rth.getADCValues().length; i++) { + double ti = HPSSVTConstants.SAMPLE_INTERVAL * i; + double fit_y = A * (Math.max(0, (ti - t0)) / constants.getTp()) * Math.exp(1 - (ti - t0) / constants.getTp()) + constants.getPedestal(); + chisq += Math.pow((fit_y - rth.getADCValues()[i]) / constants.getNoise(), 2); + } + fit.setChiSq(chisq); + + if (A > 0) { + return chisq / (rth.getADCValues().length - 2); //TODO: p-value would be better here + } else { + return Double.POSITIVE_INFINITY; + } + }
}
diff -u -r1.4 -r1.5 --- HPSShapeFitParameters.java 25 Apr 2012 18:01:32 -0000 1.4 +++ HPSShapeFitParameters.java 28 Apr 2012 19:35:09 -0000 1.5 @@ -12,127 +12,133 @@
*/ public class HPSShapeFitParameters implements GenericObject {
- private double _t0 = Double.NaN; - private double _t0Err = Double.NaN; - private double _amp = Double.NaN; - private double _ampErr = Double.NaN; - private double _tp = Double.NaN; - private double _tpErr = Double.NaN; - private double _chiSq=Double.NaN; - - public HPSShapeFitParameters() { - } - - public HPSShapeFitParameters(double t0, double amplitude) { - _t0 = t0; - _amp = amplitude; - } - - public void setT0(double t0) { - _t0 = t0; - } - - public void setAmp(double amp) { - _amp = amp; - } - - 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 void setChiSq(double _chiSq) { - this._chiSq = _chiSq; - } - - public double getT0() { - return _t0; - } - - public double getAmp() { - return _amp; - } - - public double getTp() { - return _tp; - } - - public double getT0Err() { - return _t0Err; - } - - public double getAmpErr() { - return _ampErr; - } - - public double getTpErr() { - return _tpErr; - } - public double getChiSq() { - return _chiSq; - } - - @Override - public int getNInt() { - return 0; - } - - @Override - public int getNFloat() { - return 0; - } - - @Override - public int getNDouble() { - return 7; - } - - @Override - public int getIntVal(int index) { - throw new UnsupportedOperationException("No int values in " + this.getClass().getSimpleName()); - } - - @Override - public float getFloatVal(int index) { - throw new UnsupportedOperationException("No int values in " + this.getClass().getSimpleName()); - } - - @Override - public double getDoubleVal(int index) { - switch (index) { - case 0: - return _t0; - case 1: - return _t0Err; - case 2: - return _amp; - case 3: - return _ampErr; - case 4: - return _tp; - case 5: - return _tpErr; - case 6: - return _chiSq; - default: - throw new UnsupportedOperationException("Only 7 double values in " + this.getClass().getSimpleName()); - } - - } - - @Override - public boolean isFixedSize() { - return true; - }
+ private double _t0 = Double.NaN; + private double _t0Err = Double.NaN; + private double _amp = Double.NaN; + private double _ampErr = Double.NaN; + private double _tp = Double.NaN; + private double _tpErr = Double.NaN; + private double _chiSq = Double.NaN; + + public HPSShapeFitParameters() { + } + + public HPSShapeFitParameters(double t0, double amplitude) { + _t0 = t0; + _amp = amplitude; + } + + public void setT0(double t0) { + _t0 = t0; + } + + public void setAmp(double amp) { + _amp = amp; + } + + 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 void setChiSq(double _chiSq) { + this._chiSq = _chiSq; + } + + public double getT0() { + return _t0; + } + + public double getAmp() { + return _amp; + } + + public double getTp() { + return _tp; + } + + public double getT0Err() { + return _t0Err; + } + + public double getAmpErr() { + return _ampErr; + } + + public double getTpErr() { + return _tpErr; + } + + public double getChiSq() { + return _chiSq; + } + + @Override + public int getNInt() { + return 0; + } + + @Override + public int getNFloat() { + return 0; + } + + @Override + public int getNDouble() { + return 7; + } + + @Override + public int getIntVal(int index) { + throw new UnsupportedOperationException("No int values in " + this.getClass().getSimpleName()); + } + + @Override + public float getFloatVal(int index) { + throw new UnsupportedOperationException("No int values in " + this.getClass().getSimpleName()); + } + + @Override + public double getDoubleVal(int index) { + switch (index) { + case 0: + return _t0; + case 1: + return _t0Err; + case 2: + return _amp; + case 3: + return _ampErr; + case 4: + return _tp; + case 5: + return _tpErr; + case 6: + return _chiSq; + default: + throw new UnsupportedOperationException("Only 7 double values in " + this.getClass().getSimpleName()); + } + + } + + @Override + public boolean isFixedSize() { + return true; + } + + @Override + public String toString() { + return String.format("chisq=%f\tA=%f\tAerr=%f\tT0=%f\tT0err=%f", _chiSq, _amp, _ampErr, _t0, _t0Err); + }
}
diff -u -r1.4 -r1.5 --- HPSSVTCalibrationConstants.java 25 Apr 2012 18:01:32 -0000 1.4 +++ HPSSVTCalibrationConstants.java 28 Apr 2012 19:35:09 -0000 1.5 @@ -101,17 +101,27 @@
} }
- public static ChannelConstants getChannelConstants(SiSensor sensor, int channel) { - Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel);
+ public static ChannelConstants getChannelConstants(SiSensor sensor, int channel) { + Pair<SiSensor, Integer> sensorChannel = new Pair(sensor, channel);
ChannelConstants constants = new ChannelConstants(); // ...don't have a constants file yet! noise.put(sensorChannel,20.0); pedestal.put(sensorChannel,1638.0);
- tShaping.put(sensorChannel,36.0);
+ tShaping.put(sensorChannel,35.0);
////
- constants.setNoise(noise.get(sensorChannel)); - constants.setPedestal(pedestal.get(sensorChannel)); - constants.setTp(tShaping.get(sensorChannel));
+ Double value; + value = noise.get(sensorChannel); + if (value != null) { + constants.setNoise(value); + } + value = pedestal.get(sensorChannel); + if (value != null) { + constants.setPedestal(value); + } + value = tShaping.get(sensorChannel); + if (value != null) { + constants.setTp(value); + }
return constants; }
diff -u -r1.15 -r1.16 --- TrackerReconDriver.java 17 Apr 2012 18:13:42 -0000 1.15 +++ TrackerReconDriver.java 28 Apr 2012 19:35:09 -0000 1.16 @@ -23,7 +23,7 @@
* with the {@link TrackerDigiDriver} digitization Driver. * * @author jeremym
- * @version $Id: TrackerReconDriver.java,v 1.15 2012/04/17 18:13:42 jeremy Exp $
+ * @version $Id: TrackerReconDriver.java,v 1.16 2012/04/28 19:35:09 meeg Exp $
*/ public final class TrackerReconDriver extends Driver {
@@ -195,6 +195,7 @@
hthdriver.setMaxSeperation(stripMaxSeparation); hthdriver.setTolerance(stripTolerance); // user parameter? hthdriver.setTransformToTracking(true);
+ hthdriver.setDebug(true);
add(hthdriver); //
diff -u -r1.6 -r1.7 --- HPSRawTrackerHitFitterDriver.java 25 Apr 2012 18:01:32 -0000 1.6 +++ HPSRawTrackerHitFitterDriver.java 28 Apr 2012 19:35:09 -0000 1.7 @@ -15,59 +15,61 @@
*/ public class HPSRawTrackerHitFitterDriver extends Driver {
- private HPSShaperFitAlgorithm _shaper = new DumbShaperFit(); - private String rawHitCollectionName = "SVTRawTrackerHits"; - private String fitCollectionName = "SVTShapeFitParameters"; - private String fittedHitCollectionName = "SVTFittedRawTrackerHits"; - private int genericObjectFlags = 1 << LCIOConstants.GOBIT_FIXED; - private int relationFlags = 0; - - public void setFitAlgorithm(String fitAlgorithm) { - if (fitAlgorithm.equals("Analytic")) { - _shaper = new HPSShaperAnalyticFitAlgorithm(); - } - } - - public void setFitCollectionName(String fitCollectionName) { - this.fitCollectionName = fitCollectionName; - } - - public void setFittedHitCollectionName(String fittedHitCollectionName) { - this.fittedHitCollectionName = fittedHitCollectionName; - } - - public void setRawHitCollectionName(String rawHitCollectionName) { - this.rawHitCollectionName = rawHitCollectionName; - } - - @Override - public void startOfData() { - if (rawHitCollectionName == null) { - throw new RuntimeException("The parameter ecalCollectionName was not set!"); - } - HPSSVTCalibrationConstants calibConstants=new HPSSVTCalibrationConstants(); - calibConstants.loadCalibrationConstants("Foobar"); - } - - @Override - public void process(EventHeader event) { - List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, rawHitCollectionName); - if (rawHits == null) { - throw new RuntimeException("Event is missing Raw hits collection!"); - } - List<HPSFittedRawTrackerHit> hits = new ArrayList<HPSFittedRawTrackerHit>(); - List<HPSShapeFitParameters> fits = new ArrayList<HPSShapeFitParameters>(); - - // Make a fitted hit from this cluster - for (RawTrackerHit hit : rawHits) { - 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); - hit.getDetectorElement().getReadout().addHit(hth); - } - event.put(fitCollectionName, fits, HPSShapeFitParameters.class, genericObjectFlags); - event.put(fittedHitCollectionName, hits, HPSFittedRawTrackerHit.class, relationFlags);
+ private HPSShaperFitAlgorithm _shaper = new DumbShaperFit(); + private String rawHitCollectionName = "SVTRawTrackerHits"; + private String fitCollectionName = "SVTShapeFitParameters"; + private String fittedHitCollectionName = "SVTFittedRawTrackerHits"; + private int genericObjectFlags = 1 << LCIOConstants.GOBIT_FIXED; + private int relationFlags = 0; + + public void setFitAlgorithm(String fitAlgorithm) { + if (fitAlgorithm.equals("Analytic")) { + _shaper = new HPSShaperAnalyticFitAlgorithm(); + } else { + throw new RuntimeException("Unrecognized fitAlgorithm: " + fitAlgorithm);
}
+ } + + public void setFitCollectionName(String fitCollectionName) { + this.fitCollectionName = fitCollectionName; + } + + public void setFittedHitCollectionName(String fittedHitCollectionName) { + this.fittedHitCollectionName = fittedHitCollectionName; + } + + public void setRawHitCollectionName(String rawHitCollectionName) { + this.rawHitCollectionName = rawHitCollectionName; + } + + @Override + public void startOfData() { + if (rawHitCollectionName == null) { + throw new RuntimeException("The parameter ecalCollectionName was not set!"); + } + HPSSVTCalibrationConstants calibConstants = new HPSSVTCalibrationConstants(); + calibConstants.loadCalibrationConstants("Foobar"); + } + + @Override + public void process(EventHeader event) { + List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, rawHitCollectionName); + if (rawHits == null) { + throw new RuntimeException("Event is missing SVT hits collection!"); + } + List<HPSFittedRawTrackerHit> hits = new ArrayList<HPSFittedRawTrackerHit>(); + List<HPSShapeFitParameters> fits = new ArrayList<HPSShapeFitParameters>(); + + // Make a fitted hit from this cluster + for (RawTrackerHit hit : rawHits) { + 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); + hit.getDetectorElement().getReadout().addHit(hth); + } + event.put(fitCollectionName, fits, HPSShapeFitParameters.class, genericObjectFlags); + event.put(fittedHitCollectionName, hits, HPSFittedRawTrackerHit.class, relationFlags); + }
}
diff -u -r1.4 -r1.5 --- HPSHelicalTrackHitDriver.java 19 Mar 2012 22:26:28 -0000 1.4 +++ HPSHelicalTrackHitDriver.java 28 Apr 2012 19:35:09 -0000 1.5 @@ -6,9 +6,6 @@
*/ package org.lcsim.hps.recon.tracking;
-import hep.physics.matrix.BasicMatrix; -import hep.physics.matrix.Matrix; -import hep.physics.matrix.MatrixOp;
import hep.physics.matrix.SymmetricMatrix; import hep.physics.vec.BasicHep3Vector; import hep.physics.vec.Hep3Vector;
@@ -205,7 +202,6 @@
event.put(_outname, helhits, HelicalTrackHit.class, 0); event.put(_hitrelname, hitrelations, LCRelation.class, 0); event.put(_mcrelname, mcrelations, LCRelation.class, 0);
- return;
} /**
@@ -218,32 +214,26 @@
public void OutputCollection(String outname) { _outname = outname;
- return;
} public void setOutputCollectionName(String outname) { _outname = outname;
- return;
} public void HitRelationName(String hitrelname) { _hitrelname = hitrelname;
- return;
} public void MCRelationName(String mcrelname) { _mcrelname = mcrelname;
- return;
} public void setMaxSeperation(double maxsep) { _crosser.setMaxSeparation(maxsep);
- return;
} public void setTolerance(double tolerance) { _crosser.setTolerance(tolerance);
- return;
} public void setStereoPair(String detname, int lyr1, int lyr2) {
@@ -252,7 +242,6 @@
public void setTransformToTracking(boolean xfrm) { _doTranformToTracking = xfrm;
- return;
} private List<MCParticle> getMCParticles(TrackerCluster cluster) {
@@ -404,4 +393,8 @@
private SymmetricMatrix getCovTrkSystem(SymmetricMatrix cov) { return _detToTrk.transformCovarianceToTracking(cov); }
+ + public void setDebug(boolean debug) { + this._debug = debug; + }
}
diff -u -r1.2 -r1.3 --- OnlineTracking.lcsim 28 Apr 2012 15:25:41 -0000 1.2 +++ OnlineTracking.lcsim 28 Apr 2012 19:35:09 -0000 1.3 @@ -18,8 +18,11 @@
</driver> <driver name="HPSSVTDAQMaps" type="org.lcsim.hps.recon.tracking.HPSSVTDAQMaps"/> <driver name="SVTDataToRawTrackerHitDriver" type="org.lcsim.hps.recon.tracking.SVTDataToRawTrackerHitDriver"/>
- <driver name="HPSRawTrackerHitFitterDriver" type="org.lcsim.hps.recon.tracking.HPSRawTrackerHitFitterDriver"/>
+ <driver name="HPSRawTrackerHitFitterDriver" type="org.lcsim.hps.recon.tracking.HPSRawTrackerHitFitterDriver"> + <fitAlgorithm>Analytic</fitAlgorithm> + </driver>
<driver name="TrackerHitDriver" type="org.lcsim.hps.recon.tracking.TrackerHitDriver">
+ <debug>true</debug>
<rawTrackerHitCollectionName>SVTRawTrackerHits</rawTrackerHitCollectionName> </driver> <driver name="TrackerReconDriver"
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