Author: [log in to unmask]
Date: Tue Apr 21 12:41:40 2015
New Revision: 2765
Log:
use configurable pulse shape; using dummy pulse shape parameters until HPSJAVA-495 is resolved
Added:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/PulseShape.java
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/DumbShaperFit.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperFitAlgorithm.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/DumbShaperFit.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/DumbShaperFit.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/DumbShaperFit.java Tue Apr 21 12:41:40 2015
@@ -26,7 +26,7 @@
}
@Override
- public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth) {
+ public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, PulseShape shape) {
short[] samples = rth.getADCValues();
HpsSiSensor sensor =(HpsSiSensor) rth.getDetectorElement();
int channel = rth.getIdentifierFieldValue("strip");
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java Tue Apr 21 12:41:40 2015
@@ -178,6 +178,14 @@
cluster_seeds.add(channel_number);
}
}
+
+// System.out.println("Hits to be clustered:");
+// for (int i = 0; i < 639; i++) {
+// FittedRawTrackerHit base_hit = channel_to_hit.get(i);
+// if (base_hit != null) {
+// System.out.format("channel %d\tt0 %f\t amp %f\n", i, base_hit.getT0(), base_hit.getAmp());
+// }
+// }
// Create a list of clusters
List<List<FittedRawTrackerHit>> cluster_list = new ArrayList<List<FittedRawTrackerHit>>();
@@ -259,6 +267,12 @@
} // End of loop over seeds
+// for (List<FittedRawTrackerHit> cluster : cluster_list) {
+// System.out.format("cluster with %d hits\n",cluster.size());
+// for (FittedRawTrackerHit base_hit : cluster) {
+// System.out.format("channel %d\tt0 %f\t amp %f\n", base_hit.getRawTrackerHit().getIdentifierFieldValue("strip"), base_hit.getT0(), base_hit.getAmp());
+// }
+// }
// Finished finding clusters
return cluster_list;
}
Added: java/trunk/tracking/src/main/java/org/hps/recon/tracking/PulseShape.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/PulseShape.java (added)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/PulseShape.java Tue Apr 21 12:41:40 2015
@@ -0,0 +1,78 @@
+package org.hps.recon.tracking;
+
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+
+/**
+ *
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: $
+ */
+public abstract class PulseShape {
+
+ public abstract void setParameters(int channel, HpsSiSensor sensor);
+
+ public abstract double getAmplitudePeakNorm(double time);
+
+ public abstract double getAmplitudeIntegralNorm(double time);
+
+// public abstract double getAmplitude(double time, int channel, HpsSiSensor sensor);
+ public static class CRRC extends PulseShape {
+
+ private double tp;
+
+ @Override
+ public void setParameters(int channel, HpsSiSensor sensor) {
+ tp = 50.0;
+// tp = sensor.getShapeFitParameters(channel)[HpsSiSensor.TP_INDEX];
+ }
+
+ @Override
+ public double getAmplitudePeakNorm(double time) {
+ if (time < 0) {
+ return 0;
+ }
+ return (time / tp) * Math.exp(1 - time / tp);
+ }
+
+ @Override
+ public double getAmplitudeIntegralNorm(double time) {
+ if (time < 0) {
+ return 0;
+ }
+ return (time / Math.pow(tp, 2)) * Math.exp(-time / tp);
+ }
+ }
+
+ public static class FourPole extends PulseShape {
+
+ private double tp;
+ private double tp2;
+ private double peak_t, peak_amp;
+
+ @Override
+ public void setParameters(int channel, HpsSiSensor sensor) {
+// tp = sensor.getShapeFitParameters(channel)[HpsSiSensor.TP_INDEX];
+ tp = 80.0;
+ tp2 = 12.0;
+ peak_t = 3.0 * Math.pow(tp * Math.pow(tp2, 3), 0.25); //approximate solution to exp(x)=1+x+x^2*tp/(2*tp2), where x=(1/tp2-1/tp)*t
+ peak_amp = getAmplitudeIntegralNorm(peak_t);
+ }
+
+ @Override
+ public double getAmplitudeIntegralNorm(double time) {
+ if (time < 0) {
+ return 0;
+ }
+ //===> return (time / channelConstants.getTp()) * Math.exp(1 - time / channelConstants.getTp());
+ return (Math.pow(tp, 2) / Math.pow(tp - tp2, 3))
+ * (Math.exp(-time / tp)
+ - Math.exp(-time / tp2) * (1 + time * (tp - tp2) / (tp * tp2) + 0.5 * Math.pow(time * (tp - tp2) / (tp * tp2), 2)));
+// return (time / tp) * Math.exp(1 - time / tp);
+ }
+
+ @Override
+ public double getAmplitudePeakNorm(double time) {
+ return getAmplitudeIntegralNorm(time) / peak_amp;
+ }
+ }
+}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/RawTrackerHitFitterDriver.java Tue Apr 21 12:41:40 2015
@@ -2,19 +2,14 @@
import java.util.ArrayList;
import java.util.List;
-import java.util.logging.Level;
-import java.util.logging.Logger;
//===> import org.hps.conditions.deprecated.HPSSVTCalibrationConstants.ChannelConstants;
-
-
import org.hps.readout.ecal.ReadoutTimestamp;
import org.hps.readout.svt.HPSSVTConstants;
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
//===> 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;
@@ -27,7 +22,8 @@
public class RawTrackerHitFitterDriver extends Driver {
private boolean debug = false;
- private ShaperFitAlgorithm _shaper = new DumbShaperFit();
+ private ShaperFitAlgorithm fitter = new DumbShaperFit();
+ private PulseShape shape = new PulseShape.FourPole();
private String rawHitCollectionName = "SVTRawTrackerHits";
private String fitCollectionName = "SVTShapeFitParameters";
private String fittedHitCollectionName = "SVTFittedRawTrackerHits";
@@ -65,15 +61,25 @@
public void setFitAlgorithm(String fitAlgorithm) {
if (fitAlgorithm.equals("Analytic")) {
- _shaper = new ShaperAnalyticFitAlgorithm();
+ fitter = new ShaperAnalyticFitAlgorithm();
} else if (fitAlgorithm.equals("Linear")) {
- _shaper = new ShaperLinearFitAlgorithm(1);
+ fitter = new ShaperLinearFitAlgorithm(1);
} else if (fitAlgorithm.equals("PileupAlways")) {
- _shaper = new ShaperPileupFitAlgorithm(1.0);
+ fitter = new ShaperPileupFitAlgorithm(1.0);
} else if (fitAlgorithm.equals("Pileup")) {
- _shaper = new ShaperPileupFitAlgorithm();
+ fitter = new ShaperPileupFitAlgorithm();
} else {
throw new RuntimeException("Unrecognized fitAlgorithm: " + fitAlgorithm);
+ }
+ }
+
+ public void setPulseShape(String pulseShape) {
+ if (pulseShape.equals("CR-RC")) {
+ shape = new PulseShape.CRRC();
+ } else if (pulseShape.equals("FourPole")) {
+ shape = new PulseShape.FourPole();
+ } else {
+ throw new RuntimeException("Unrecognized pulseShape: " + pulseShape);
}
}
@@ -91,7 +97,7 @@
@Override
public void startOfData() {
- _shaper.setDebug(debug);
+ fitter.setDebug(debug);
if (rawHitCollectionName == null) {
throw new RuntimeException("The parameter ecalCollectionName was not set!");
}
@@ -117,7 +123,7 @@
HpsSiSensor sensor = (HpsSiSensor) hit.getDetectorElement();
//===> ChannelConstants constants = HPSSVTCalibrationConstants.getChannelConstants((SiSensor) hit.getDetectorElement(), strip);
//for (ShapeFitParameters fit : _shaper.fitShape(hit, constants)) {
- for (ShapeFitParameters fit : _shaper.fitShape(hit)) {
+ for (ShapeFitParameters fit : fitter.fitShape(hit, shape)) {
if (correctT0Shift) {
//===> fit.setT0(fit.getT0() - constants.getT0Shift());
fit.setT0(fit.getT0() - sensor.getT0Shift());
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java Tue Apr 21 12:41:40 2015
@@ -24,7 +24,7 @@
}
@Override
- public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth) {
+ public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, PulseShape shape) {
short[] samples = rth.getADCValues();
HpsSiSensor sensor =(HpsSiSensor) rth.getDetectorElement();
int channel = rth.getIdentifierFieldValue("strip");
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperFitAlgorithm.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperFitAlgorithm.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperFitAlgorithm.java Tue Apr 21 12:41:40 2015
@@ -12,7 +12,7 @@
public interface ShaperFitAlgorithm {
//===> public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, ChannelConstants constants);
- public Collection<ShapeFitParameters> fitShape(RawTrackerHit rawHit);
+ public Collection<ShapeFitParameters> fitShape(RawTrackerHit rawHit, PulseShape shape);
public void setDebug(boolean debug);
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java Tue Apr 21 12:41:40 2015
@@ -1,11 +1,8 @@
package org.hps.recon.tracking;
-import java.io.OutputStream;
-import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collection;
import java.util.logging.Level;
-import java.util.logging.LogManager;
import java.util.logging.Logger;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
@@ -38,14 +35,15 @@
final double[] amplitudes;
final double[] amplitudeErrors;
//===> private ChannelConstants channelConstants;
- HpsSiSensor sensor = null;
+ private HpsSiSensor sensor;
+ private int channel;
+ private PulseShape shape;
private final double[] sigma = new double[HPSSVTConstants.TOTAL_NUMBER_OF_SAMPLES];
private final double[] y = new double[HPSSVTConstants.TOTAL_NUMBER_OF_SAMPLES];
private int firstUsedSample;
private int nUsedSamples;
private int firstFittedPulse;
private int nFittedPulses;
- private int channel;
private boolean debug = false;
private static final Logger minuitLoggger = Logger.getLogger("org.freehep.math.minuit");
@@ -67,19 +65,19 @@
@Override
//===> public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, ChannelConstants constants) {
- public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth) {
+ public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, PulseShape shape) {
short[] samples = rth.getADCValues();
- HpsSiSensor sensor =(HpsSiSensor) rth.getDetectorElement();
- int channel = rth.getIdentifierFieldValue("strip");
- return fitShape(channel, samples, sensor);
+ sensor = (HpsSiSensor) rth.getDetectorElement();
+ channel = rth.getIdentifierFieldValue("strip");
+ this.shape = shape;
+ shape.setParameters(channel, sensor);
+ return fitShape(samples);
//===> return this.fitShape(rth.getADCValues(), constants);
}
- public Collection<ShapeFitParameters> fitShape(int channelNumber, short[] samples, HpsSiSensor siSensor) {
- //===> public Collection<ShapeFitParameters> fitShape(short[] samples, ChannelConstants constants) {
+ public Collection<ShapeFitParameters> fitShape(short[] samples) {
+ //===> public Collection<ShapeFitParameters> fitShape(short[] samples, ChannelConstants constants) {
// channelConstants = constants;
- this.sensor = siSensor;
- this.channel = channelNumber;
double[] signal = new double[HPSSVTConstants.TOTAL_NUMBER_OF_SAMPLES];
for (int i = 0; i < samples.length; i++) {
@@ -127,14 +125,30 @@
fit.setT0(min.userState().value(i));
fit.setT0Err(min.userState().error(i));
-// if (min.isValid()) {
-// MnMinos minos = new MnMinos(this, min);
-// MinosError t0err = minos.minos(0);
-// if (t0err.isValid()) {
+
+// MinosError t0err = null;
+// if (min.isValid() && min.edm() > 0) {
+// MnMinos minos = null;
+//
+// try {
+// minos = new MnMinos(this, min);
+// t0err = minos.minos(0);
+// } catch (RuntimeException e) {
+// if (debug) {
+// System.out.println(e);
+// }
+// }
+// }
+// if (t0err != null && t0err.isValid()) {
+// if (debug) {
+// System.out.format("fitter error %f, minos lower %f, upper %f\n", min.userState().error(i), t0err.lower(), t0err.upper());
+// }
// fit.setT0Err((t0err.lower() + t0err.upper()) / 2);
+// } else {
+// if (debug) {
+// System.out.format("fitter error %f\n", min.userState().error(i));
+// }
// }
-// }
-
// System.out.println(fit);
fits.add(fit);
}
@@ -171,7 +185,7 @@
//subtract first pulse from fit input
for (int i = 0; i < samples.length; i++) {
//===> fitData[i] -= amplitudes[firstFittedPulse] * getAmplitude(HPSSVTConstants.SAMPLING_INTERVAL * i - frontFit.userState().value(0), channelConstants);
- fitData[i] -= amplitudes[firstFittedPulse] * getAmplitude(HPSSVTConstants.SAMPLING_INTERVAL * i - frontFit.userState().value(0), this.channel, this.sensor);
+ fitData[i] -= amplitudes[firstFittedPulse] * shape.getAmplitudePeakNorm(HPSSVTConstants.SAMPLING_INTERVAL * i - frontFit.userState().value(0));
}
if (debug) {
@@ -315,7 +329,7 @@
for (int j = 0; j < nUsedSamples; j++) {
for (int i = 0; i < nFittedPulses; i++) {
//===> sc_mat.setEntry(i, j, getAmplitude(HPSSVTConstants.SAMPLING_INTERVAL * (firstUsedSample + j) - times[i], channelConstants) / sigma[firstUsedSample + j]);
- sc_mat.setEntry(i, j, getAmplitude(HPSSVTConstants.SAMPLING_INTERVAL * (firstUsedSample + j) - times[i], this.channel, this.sensor) / sigma[firstUsedSample + j]);
+ sc_mat.setEntry(i, j, shape.getAmplitudePeakNorm(HPSSVTConstants.SAMPLING_INTERVAL * (firstUsedSample + j) - times[i]) / sigma[firstUsedSample + j]);
}
y_vec.setEntry(j, y[firstUsedSample + j] / sigma[firstUsedSample + j]);
var_vec.setEntry(j, sigma[firstUsedSample + j] * sigma[firstUsedSample + j]);
@@ -351,16 +365,15 @@
return chisq;
}
- //===> private static double getAmplitude(double time, ChannelConstants channelConstants) {
- private static double getAmplitude(double time, int channel, HpsSiSensor sensor) {
- if (time < 0) {
- return 0;
- }
- double tp = sensor.getShapeFitParameters(channel)[HpsSiSensor.TP_INDEX];
- //===> return (time / channelConstants.getTp()) * Math.exp(1 - time / channelConstants.getTp());
- return (time / tp) * Math.exp(1 - time / tp);
- }
-
+// //===> private static double getAmplitude(double time, ChannelConstants channelConstants) {
+// private static double getAmplitude(double time, int channel, HpsSiSensor sensor) {
+// if (time < 0) {
+// return 0;
+// }
+// double tp = sensor.getShapeFitParameters(channel)[HpsSiSensor.TP_INDEX];
+// //===> return (time / channelConstants.getTp()) * Math.exp(1 - time / channelConstants.getTp());
+// return (time / tp) * Math.exp(1 - time / tp);
+// }
@Override
public double valueOf(double[] times) {
return doLinFit(times);
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperPileupFitAlgorithm.java Tue Apr 21 12:41:40 2015
@@ -26,15 +26,15 @@
}
//===> public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, HPSSVTCalibrationConstants.ChannelConstants constants) {
- public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth) {
+ public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, PulseShape shape) {
//===> Collection<ShapeFitParameters> fittedPulses = onePulseFitter.fitShape(rth, constants);
- Collection<ShapeFitParameters> fittedPulses = onePulseFitter.fitShape(rth);
+ Collection<ShapeFitParameters> fittedPulses = onePulseFitter.fitShape(rth, shape);
double singlePulseChiProb = fittedPulses.iterator().next().getChiProb();
totalFits++;
if (singlePulseChiProb < refitThreshold) {
refitAttempts++;
//===> Collection<ShapeFitParameters> doublePulse = twoPulseFitter.fitShape(rth, constants);
- Collection<ShapeFitParameters> doublePulse = twoPulseFitter.fitShape(rth);
+ Collection<ShapeFitParameters> doublePulse = twoPulseFitter.fitShape(rth, shape);
double doublePulseChiProb = doublePulse.iterator().next().getChiProb();
if (doublePulseChiProb > singlePulseChiProb) {
refitsAccepted++;
|
