8 modified files
java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt
--- java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt/SVTMonitoringPlots.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt/SVTMonitoringPlots.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -6,6 +6,7 @@
import hep.aida.IProfile1D;
import java.util.List;
+import org.apache.commons.math3.special.Gamma;
import org.hps.conditions.deprecated.HPSSVTCalibrationConstants;
import org.hps.conditions.deprecated.SvtUtils;
@@ -150,7 +151,7 @@
int layer = fit.getRawTrackerHit().getIdentifierFieldValue("layer"); // 1-10; axial layers are odd layers; stereo layers are even
int module = fit.getRawTrackerHit().getIdentifierFieldValue("module"); // 0-1; module number is top or bottom
int strip = fit.getRawTrackerHit().getIdentifierFieldValue("strip");
- if (fit.getShapeFitParameters().getChiSq() < 5) {
+ if (fit.getShapeFitParameters().getChiProb() > Gamma.regularizedGammaQ(4, 5)) {
double noise = HPSSVTCalibrationConstants.getNoise((SiSensor) fit.getRawTrackerHit().getDetectorElement(), strip);
if (fit.getAmp() > 4 * noise) {
t0s[module][layer - 1].fill(strip, fit.getT0());
@@ -162,6 +163,6 @@
}
public void endOfData() {
- //plotterFrame.dispose();
+ //plotterFrame.dispose();
}
-}
\ No newline at end of file
+}
java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt
--- java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt/SVTPulseFitPlots.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/monitoring-drivers/src/main/java/org/hps/monitoring/drivers/svt/SVTPulseFitPlots.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -42,7 +42,7 @@
private String outputPlots = null;
private IHistogram1D[][] t0 = new IHistogram1D[2][10];
private IHistogram1D[][] amp = new IHistogram1D[2][10];
- private IHistogram1D[][] chisq = new IHistogram1D[2][10];
+ private IHistogram1D[][] chiprob = new IHistogram1D[2][10];
private IHistogram2D[][] t0a = new IHistogram2D[2][10];
private IHistogram2D[][] shape = new IHistogram2D[2][10];
// private IHistogram2D shape;
@@ -100,8 +100,8 @@
plotter.region(region).plot(t0[module][layer]);
amp[module][layer] = aida.histogram1D(sensor.getName() + "_amplitude", 50, 0, 2000.0);
plotter2.region(region).plot(amp[module][layer]);
- chisq[module][layer] = aida.histogram1D(sensor.getName() + "_chisq", 100, 0, 10.0);
- plotter3.region(region).plot(chisq[module][layer]);
+ chiprob[module][layer] = aida.histogram1D(sensor.getName() + "_chiprob", 100, 0, 1.0);
+ plotter3.region(region).plot(chiprob[module][layer]);
t0a[module][layer] = aida.histogram2D(sensor.getName() + " A vs. T0", 100, -100, 100, 100, 0, 2000);
plotter4.region(region).plot(t0a[module][layer]);
shape[module][layer] = aida.histogram2D(sensor.getName() + " Shape", 200, -1, 3, 200, -0.5, 2);
@@ -136,7 +136,7 @@
String sensorName = sensor.getName();
aida.histogram1D(sensorName + "_timing").fill(fittedT0);
aida.histogram1D(sensorName + "_amplitude").fill(fittedAmp);
- aida.histogram1D(sensorName + "_chisq").fill(fit.getShapeFitParameters().getChiSq());
+ aida.histogram1D(sensorName + "_chiprob").fill(fit.getShapeFitParameters().getChiProb());
double noise = HPSSVTCalibrationConstants.getNoise(sensor, strip);
double pedestal = HPSSVTCalibrationConstants.getPedestal(sensor, strip);
@@ -167,7 +167,7 @@
for (int layer = 0; layer < 10; layer++) {
t0[module][layer].reset();
amp[module][layer].reset();
- chisq[module][layer].reset();
+ chiprob[module][layer].reset();
t0a[module][layer].reset();
shape[module][layer].reset();
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/NearestNeighborRMSClusterer.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -7,7 +7,7 @@
import java.util.List;
import java.util.Map;
import java.util.Set;
-
+import org.apache.commons.math3.special.Gamma;
import org.hps.conditions.deprecated.HPSSVTCalibrationConstants;
import org.hps.conditions.deprecated.HPSSVTConstants;
import org.lcsim.detector.identifier.IIdentifier;
@@ -28,7 +28,7 @@
private double _cluster_threshold;
private double _meanTime = 24;
private double _timeWindow = 48;
- private double _maxChisq = 20.0;
+ private final double _minChiProb = Gamma.regularizedGammaQ(4, 20);
/**
* Instantiate NearestNeighborRMS with specified thresholds. Seed threshold is the minimum
@@ -231,7 +231,7 @@
}
private boolean passChisqCut(FittedRawTrackerHit hit) {
- return hit.getShapeFitParameters().getChiSq() < _maxChisq;
+ return hit.getShapeFitParameters().getChiProb() > _minChiProb;
}
public int getNeighborCell(int cell, int ncells_0, int ncells_1) {
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShapeFitParameters.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShapeFitParameters.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -14,7 +14,7 @@
private double _ampErr = Double.NaN;
// private double _tp = Double.NaN;
// private double _tpErr = Double.NaN;
- private double _chiSq = Double.NaN;
+ private double _chiProb = Double.NaN;
public ShapeFitParameters() {
}
@@ -48,8 +48,8 @@
// this._tpErr = _tpErr;
// }
- public void setChiSq(double _chiSq) {
- this._chiSq = _chiSq;
+ public void setChiProb(double _chiProb) {
+ this._chiProb = _chiProb;
}
public double getT0() {
@@ -76,8 +76,8 @@
// return _tpErr;
// }
- public double getChiSq() {
- return _chiSq;
+ public double getChiProb() {
+ return _chiProb;
}
public static double getT0(GenericObject object) {
@@ -104,7 +104,7 @@
// return object.getDoubleVal(5);
// }
- public static double getChisq(GenericObject object) {
+ public static double getChiProb(GenericObject object) {
return object.getDoubleVal(6);
}
@@ -149,9 +149,9 @@
// case 5:
// return _tpErr;
case 4:
- return _chiSq;
+ return _chiProb;
default:
- throw new UnsupportedOperationException("Only 7 double values in " + this.getClass().getSimpleName());
+ throw new UnsupportedOperationException("Only 5 double values in " + this.getClass().getSimpleName());
}
}
@@ -163,6 +163,6 @@
@Override
public String toString() {
- return String.format("chisq=%f\tA=%f\tAerr=%f\tT0=%f\tT0err=%f", _chiSq, _amp, _ampErr, _t0, _t0Err);
+ return String.format("chiprob=%f\tA=%f\tAerr=%f\tT0=%f\tT0err=%f", _chiProb, _amp, _ampErr, _t0, _t0Err);
}
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperAnalyticFitAlgorithm.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -2,23 +2,22 @@
import java.util.ArrayList;
import java.util.Collection;
+import org.apache.commons.math3.special.Gamma;
import org.hps.conditions.deprecated.HPSSVTCalibrationConstants.ChannelConstants;
import org.hps.conditions.deprecated.HPSSVTConstants;
import org.lcsim.event.RawTrackerHit;
-//import org.lcsim.math.chisq.ChisqProb;
/**
- * 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 Sho Uemura
*/
public class ShaperAnalyticFitAlgorithm implements ShaperFitAlgorithm {
@Override
public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, ChannelConstants constants) {
- short[] samples = rth.getADCValues();
- return this.fitShape(samples, constants);
+ return this.fitShape(rth.getADCValues(), constants);
}
public Collection<ShapeFitParameters> fitShape(short[] samples, ChannelConstants constants) {
@@ -98,10 +97,10 @@
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 - samples[i]) / constants.getNoise(), 2);
}
- fit.setChiSq(chisq);
- if (A > 0) {
+ if (A > 0 && chisq < Double.POSITIVE_INFINITY) {
// return ChisqProb.gammp(samples.length - 2, chisq);
+ fit.setChiProb(Gamma.regularizedGammaQ(samples.length - 2, chisq));
return chisq / (samples.length - 2);
} else {
return Double.POSITIVE_INFINITY;
java/trunk/tracking/src/main/java/org/hps/recon/tracking
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/ShaperLinearFitAlgorithm.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -5,8 +5,11 @@
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.CholeskyDecomposition;
+import org.apache.commons.math3.linear.DecompositionSolver;
+import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
+import org.apache.commons.math3.special.Gamma;
import org.freehep.math.minuit.FCNBase;
import org.freehep.math.minuit.FunctionMinimum;
import org.freehep.math.minuit.MinosError;
@@ -16,7 +19,6 @@
import org.hps.conditions.deprecated.HPSSVTCalibrationConstants.ChannelConstants;
import org.hps.conditions.deprecated.HPSSVTConstants;
import org.lcsim.event.RawTrackerHit;
-//import org.lcsim.math.chisq.ChisqProb;
/**
* Fast fitter; currently only fits single hits. Uses Tp from ChannelConstants;
@@ -27,7 +29,6 @@
public class ShaperLinearFitAlgorithm implements ShaperFitAlgorithm, FCNBase {
final int nPeaks;
- final double[] times;
final double[] amplitudes;
final double[] amplitudeErrors;
private ChannelConstants channelConstants;
@@ -37,22 +38,19 @@
public ShaperLinearFitAlgorithm() {
nPeaks = 1;
- times = new double[nPeaks];
amplitudes = new double[nPeaks];
amplitudeErrors = new double[nPeaks];
}
public ShaperLinearFitAlgorithm(int nPeaks) {
this.nPeaks = nPeaks;
- times = new double[nPeaks];
amplitudes = new double[nPeaks];
amplitudeErrors = new double[nPeaks];
}
@Override
public Collection<ShapeFitParameters> fitShape(RawTrackerHit rth, ChannelConstants constants) {
- short[] samples = rth.getADCValues();
- return this.fitShape(samples, constants);
+ return this.fitShape(rth.getADCValues(), constants);
}
public Collection<ShapeFitParameters> fitShape(short[] samples, ChannelConstants constants) {
@@ -61,7 +59,7 @@
sigma = new double[samples.length];
usedSamples = new int[samples.length];
for (int i = 0; i < samples.length; i++) {
- y[i] = samples[i];
+ y[i] = samples[i] - constants.getPedestal();
sigma[i] = constants.getNoise();
usedSamples[i] = i;
}
@@ -73,22 +71,24 @@
MnSimplex simplex = new MnSimplex(this, myParams, 2);
FunctionMinimum min = simplex.minimize();
- MnMinos minos = new MnMinos(this, min);
- MinosError t0err = minos.minos(0);
-
ShapeFitParameters fit = new ShapeFitParameters();
fit.setAmp(amplitudes[0]);
- fit.setAmpErr(amplitudes[0]);
- fit.setChiSq(min.fval());
- fit.setT0(times[0]);
- fit.setT0Err((t0err.lower() + t0err.upper()) / 2);
+ fit.setAmpErr(amplitudeErrors[0]);
+ fit.setChiProb(Gamma.regularizedGammaQ(samples.length - 2, min.fval()));
+
+ fit.setT0(min.userState().value(0));
- // System.out.format("%f\t%f\t%f\t%f\t%f\t%f\n", samples[0] - constants.getPedestal(),
- // samples[1] - constants.getPedestal(), samples[2] - constants.getPedestal(), samples[3] -
- // constants.getPedestal(), samples[4] - constants.getPedestal(), samples[5] -
- // constants.getPedestal());
- // System.out.println("start = " + bestStart + ", " + fit);
+ fit.setT0Err(HPSSVTConstants.SAMPLING_INTERVAL);
+// if (min.isValid()) {
+// MnMinos minos = new MnMinos(this, min);
+// MinosError t0err = minos.minos(0);
+// if (t0err.isValid()) {
+// fit.setT0Err((t0err.lower() + t0err.upper()) / 2);
+// }
+// }
+
+// System.out.println(fit);
ArrayList<ShapeFitParameters> fits = new ArrayList<ShapeFitParameters>();
fits.add(fit);
return fits;
@@ -108,14 +108,25 @@
}
RealVector a_vec = sc_mat.operate(y_vec);
RealMatrix coeff_mat = sc_mat.multiply(sc_mat.transpose());
- CholeskyDecomposition a_qr = new CholeskyDecomposition(coeff_mat);
- RealVector solved_amplitudes = a_qr.getSolver().solve(a_vec);
- RealVector amplitude_err = a_qr.getSolver().solve(var_vec);
+ DecompositionSolver a_solver;
+ RealVector solved_amplitudes, amplitude_err;
+ try {
+ CholeskyDecomposition a_cholesky = new CholeskyDecomposition(coeff_mat);
+ a_solver = a_cholesky.getSolver();
+ solved_amplitudes = a_solver.solve(a_vec);
+ amplitude_err = a_solver.solve(sc_mat.operate(var_vec));
+ } catch (NonPositiveDefiniteMatrixException e) {
+ solved_amplitudes = new ArrayRealVector(nPeaks, 0.0);
+ amplitude_err = new ArrayRealVector(nPeaks, Double.POSITIVE_INFINITY);
+ }
+
+ double chisq = y_vec.subtract(sc_mat.preMultiply(solved_amplitudes)).getNorm();
+
for (int i = 0; i < times.length; i++) {
amplitudes[i] = solved_amplitudes.getEntry(i);
amplitudeErrors[i] = Math.sqrt(amplitude_err.getEntry(i));
}
- return y_vec.subtract(sc_mat.operate(solved_amplitudes)).getNorm();
+ return chisq;
}
private static double getAmplitude(double time, ChannelConstants channelConstants) {
java/trunk/users/src/main/java/org/hps/users/omoreno
--- java/trunk/users/src/main/java/org/hps/users/omoreno/SvtPerformance.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/users/src/main/java/org/hps/users/omoreno/SvtPerformance.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -310,7 +310,6 @@
ChannelConstants constants = null;
double clusterAmplitude, maxClusterAmplitude;
double noise = 0;
- double chiSquared = -1;
double trkChiSquared = 0;
double hitTime = 0;
double hitX, hitY, pedestal;
@@ -353,7 +352,7 @@
hitsPerCluster = stripHit.rawhits().size();
noise = 0;
bad_channel = 0;
- chiSquared = -1;
+ double chiSquaredProb = -1;
for (Object rh : stripHit.rawhits()) {
RawTrackerHit rawHit = (RawTrackerHit) rh;
@@ -388,7 +387,7 @@
maxClusterAmplitude = fit.getAmp();
}
if (stripHit.rawhits().size() == 1) {
- chiSquared = fit.getChiSq();
+ chiSquaredProb = fit.getChiProb();
}
noise += Math.pow(sensor.getNoise(channel), 2);
clusterAmplitude += fit.getAmp();
@@ -407,7 +406,7 @@
performanceWriter.write(runNumber + " " + eventNumber + " 1 " + sensor.getLayerNumber() + " ");
}
performanceWriter.write(maxClusterChannel + " " + clusterAmplitude + " " + noise + " " + hitsPerCluster + " "
- + bad_channel + " " + chiSquared + " " + hitX + " " + hitY + " " + trkChiSquared + " "
+ + bad_channel + " " + chiSquaredProb + " " + hitX + " " + hitY + " " + trkChiSquared + " "
+ hitTime + "\n");
} catch (IOException exception) {
exception.printStackTrace();
java/trunk/users/src/main/java/org/hps/users/omoreno
--- java/trunk/users/src/main/java/org/hps/users/omoreno/SvtQA.java 2014-08-26 23:21:58 UTC (rev 913)
+++ java/trunk/users/src/main/java/org/hps/users/omoreno/SvtQA.java 2014-08-27 02:07:14 UTC (rev 914)
@@ -301,9 +301,9 @@
//--- Chi Squared vs Channel ---//
//------------------------------//
if (enableChiSquaredvsChannel) {
- title = sensorName + " - Chi Squared vs Channel #";
+ title = sensorName + " - Chi Squared Probability vs Channel #";
plotters.add(PlotUtils.setupPlotter(title, 0, 0));
- histo2D = aida.histogram2D(title, 640, 0, 639, 300, 0, 100);
+ histo2D = aida.histogram2D(title, 640, 0, 639, 300, 0, 1.0);
histos2D.add(histo2D);
PlotUtils.setup2DRegion(plotters.get(plotterIndex), title, 0, "Channel #", "Chi Squared", histo2D);
plotterIndex++;
@@ -558,8 +558,8 @@
// Fill Chi Squared vs Channel # plots
if (enableChiSquaredvsChannel && SvtUtils.getInstance().getDescription(sensor).equals(sensorName)) {
- title = sensorName + " - Chi Squared vs Channel #";
- aida.histogram2D(title).fill(channel, fittedHit.getShapeFitParameters().getChiSq());
+ title = sensorName + " - Chi Squared Probability vs Channel #";
+ aida.histogram2D(title).fill(channel, fittedHit.getShapeFitParameters().getChiProb());
}
}
SVNspam 0.1
