java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal
--- java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal/FADCEcalReadoutDriver.java 2014-04-21 19:03:07 UTC (rev 505)
+++ java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal/FADCEcalReadoutDriver.java 2014-04-21 19:04:13 UTC (rev 506)
@@ -43,7 +43,6 @@
private static final int ECAL_WINDOW_MODE = 1;
private static final int ECAL_PULSE_MODE = 2;
private static final int ECAL_PULSE_INTEGRAL_MODE = 3;
-
String ecalName = "Ecal";
Subdetector ecal;
//buffer for preamp signals (units of volts, no pedestal)
@@ -51,7 +50,7 @@
//ADC pipeline for readout (units of ADC counts)
private Map<Long, FADCPipeline> pipelineMap = null;
//buffer for window sums
- private Map<Long, Double> sumMap = null;
+ private Map<Long, Integer> sumMap = null;
//buffer for timestamps
private Map<Long, Integer> timeMap = null;
//queue for hits to be output to clusterer
@@ -60,10 +59,8 @@
private int bufferLength = 100;
//length of readout pipeline (in readout cycles)
private int pipelineLength = 2000;
- //switch between two pulse shape functions
- private boolean useCRRCShape = true;
- //shaper time constant in ns; negative values generate square pulses of the given width (for test run sim)
- private double tp = 14.0;
+ //shaper time constant in ns
+ private double tp = 6.95;
//delay (number of readout periods) between start of summing window and output of hit to clusterer
private int delay0 = 32;
//start of readout window relative to trigger time (in readout cycles)
@@ -86,16 +83,23 @@
//output collection name for hits read out from trigger
private String ecalReadoutCollectionName = "EcalReadoutHits";
private int mode = ECAL_PULSE_INTEGRAL_MODE;
- private int readoutThreshold = 50;
- private int triggerThreshold = 50;
- //amplitude ADC counts/GeV
-// private double gain = 0.5*1000 * 80.0 / 60;
- private double scaleFactor = 128;
+ private int readoutThreshold = 10;
+ private int triggerThreshold = 10;
+ private double scaleFactor = 1;
private double fixedGain = -1;
- private boolean constantTriggerWindow = false;
+ private boolean constantTriggerWindow = true;
private boolean addNoise = false;
private double pePerMeV = 2.0; //photoelectrons per MeV, used to calculate noise
+ //switch between test run and 2014 definitions of gain constants
+ private boolean use2014Gain = true;
+ //switch between three pulse shape functions
+ private PulseShape pulseShape = PulseShape.ThreePole;
+ public enum PulseShape {
+
+ CRRC, DoubleGaussian, ThreePole
+ }
+
public FADCEcalReadoutDriver() {
flags = 0;
flags += 1 << LCIOConstants.RCHBIT_TIME; //store timestamp
@@ -156,10 +160,14 @@
this.coincidenceWindow = coincidenceWindow;
}
- public void setUseCRRCShape(boolean useCRRCShape) {
- this.useCRRCShape = useCRRCShape;
+ public void setUse2014Gain(boolean use2014Gain) {
+ this.use2014Gain = use2014Gain;
}
+ public void setPulseShape(String pulseShape) {
+ this.pulseShape = PulseShape.valueOf(pulseShape);
+ }
+
public void setTp(double tp) {
this.tp = tp;
}
@@ -223,15 +231,17 @@
pipeline.step();
double currentValue = signalBuffer.currentValue() * ((Math.pow(2, nBit) - 1) / maxVolt); //12-bit ADC with maxVolt V range
- double pedestal = EcalConditions.physicalToPedestal(cellID);
- pipeline.writeValue(Math.min((int) Math.round(pedestal + currentValue), (int) Math.pow(2, nBit))); //ADC can't return a value larger than 4095; 4096 (overflow) is returned for any input >2V
+ int pedestal = (int) Math.round(EcalConditions.physicalToPedestal(cellID));
+ int digitizedValue = Math.min((int) Math.round(pedestal + currentValue), (int) Math.pow(2, nBit)); //ADC can't return a value larger than 4095; 4096 (overflow) is returned for any input >2V
+ pipeline.writeValue(digitizedValue);
+ int pedestalSubtractedValue = digitizedValue - pedestal;
//System.out.println(signalBuffer.currentValue() + " " + currentValue + " " + pipeline.currentValue());
- Double sum = sumMap.get(cellID);
- if (sum == null && currentValue > triggerThreshold) {
+ Integer sum = sumMap.get(cellID);
+ if (sum == null && pedestalSubtractedValue > triggerThreshold) {
timeMap.put(cellID, readoutCounter);
if (constantTriggerWindow) {
- double sumBefore = 0;
+ int sumBefore = 0;
for (int i = 0; i < numSamplesBefore; i++) {
if (debug) {
System.out.format("trigger %d, %d: %d\n", cellID, i, pipeline.getValue(numSamplesBefore - i - 1));
@@ -240,7 +250,7 @@
}
sumMap.put(cellID, sumBefore);
} else {
- sumMap.put(cellID, currentValue);
+ sumMap.put(cellID, pedestalSubtractedValue);
}
}
if (sum != null) {
@@ -259,15 +269,15 @@
sumMap.remove(cellID);
}
} else {
- if (currentValue < triggerThreshold || timeMap.get(cellID) + delay0 == readoutCounter) {
+ if (pedestalSubtractedValue < triggerThreshold || timeMap.get(cellID) + delay0 == readoutCounter) {
// System.out.printf("sum = %f\n",sum);
outputQueue.add(new HPSRawCalorimeterHit(cellID,
- (int) Math.round((sum + currentValue) / scaleFactor),
+ (int) Math.round((sum + pedestalSubtractedValue) / scaleFactor),
64 * timeMap.get(cellID),
readoutCounter - timeMap.get(cellID) + 1));
sumMap.remove(cellID);
} else {
- sumMap.put(cellID, sum + currentValue);
+ sumMap.put(cellID, sum + pedestalSubtractedValue);
}
}
}
@@ -423,8 +433,8 @@
if (addNoise) {
//add preamp noise and photoelectron Poisson noise in quadrature
double noise;
- if (!useCRRCShape) {
- noise = Math.sqrt(Math.pow(EcalConditions.physicalToNoise(hit.getCellID()) * EcalConditions.physicalToGain(hit.getCellID()) * ECalUtils.gainFactor * ECalUtils.ecalReadoutPeriod, 2) + hit.getRawEnergy() * ECalUtils.MeV / pePerMeV);
+ if (use2014Gain) {
+ noise = Math.sqrt(Math.pow(EcalConditions.physicalToNoise(hit.getCellID()) * EcalConditions.physicalToGain(hit.getCellID()) * ECalUtils.gainFactor * ECalUtils.ecalReadoutPeriod, 2) + hit.getRawEnergy() / (ECalUtils.lightYield * ECalUtils.quantumEff * ECalUtils.surfRatio));
} else {
noise = Math.sqrt(Math.pow(EcalConditions.physicalToNoise(hit.getCellID()) * EcalConditions.physicalToGain(hit.getCellID()) * ECalUtils.MeV, 2) + hit.getRawEnergy() * ECalUtils.MeV / pePerMeV);
}
@@ -439,7 +449,7 @@
@Override
protected void initReadout() {
//initialize buffers
- sumMap = new HashMap<Long, Double>();
+ sumMap = new HashMap<Long, Integer>();
timeMap = new HashMap<Long, Integer>();
outputQueue = new PriorityQueue(20, new HPSRawCalorimeterHit.TimeComparator());
resetFADCBuffers();
@@ -467,11 +477,17 @@
}
private double pulseAmplitude(double time, long cellID) {
- if (useCRRCShape) {
- if (time <= 0.0) {
- return 0.0;
+ if (use2014Gain) {
+ //if fixedGain is set, multiply the default gain by this factor
+ double corrGain;
+ if (fixedGain > 0) {
+ corrGain = fixedGain;
+ } else {
+ corrGain = 1.0 / EcalConditions.physicalToGain(cellID);
}
+ return corrGain * readoutGain * pulseAmplitude(time, pulseShape, tp);
+ } else {
//normalization constant from cal gain (MeV/integral bit) to amplitude gain (amplitude bit/GeV)
double gain;
if (fixedGain > 0) {
@@ -480,37 +496,39 @@
gain = readoutPeriod / (EcalConditions.physicalToGain(cellID) * ECalUtils.MeV * ((Math.pow(2, nBit) - 1) / maxVolt));
}
- if (tp > 0.0) {
- return gain * ((time / tp) * Math.exp(1.0 - time / tp)) / (tp * Math.E);
- } else {
- if (time < -tp) {
- return 1.0;
- } else {
- return 0.0;
- }
- }
- } else { // According to measurements the output signal can be fitted by two gaussians, one for the rise of the signal, one for the fall
- // Time corresponds to t-(t_eve+pulseDelay) such that the maximum of the amplitude is reached pulseDelay ns after the event t_eve
- // Without the coefficient, the integral is equal to 1
- if (time <= 0.0) {
- return 0.0;
- }
+ return gain * pulseAmplitude(time, pulseShape, tp);
+ }
+ }
- //if fixedGain is set, multiply the default gain by this factor
- double corrGain = 1.0;
- if (fixedGain > 0) {
- corrGain = fixedGain;
- } else {
- corrGain = 1.0 / EcalConditions.physicalToGain(cellID);
- }
+ /**
+ * Returns pulse amplitude at the given time (relative to hit time).
+ * Amplitude is normalized so the pulse integral is 1.
+ *
+ * @param time
+ * @return
+ */
+ public static double pulseAmplitude(double time, PulseShape shape, double shapingTime) {
+ if (time <= 0.0) {
+ return 0.0;
+ }
+ switch (shape) {
+ case CRRC:
+ //peak at tp
+ //peak value 1/(tp*e)
+ return ((time / (shapingTime * shapingTime)) * Math.exp(-time / shapingTime));
+ case DoubleGaussian:
+ //According to measurements the output signal can be fitted by two gaussians, one for the rise of the signal, one for the fall
+ //peak at 3*riseTime
+ //peak value 1/norm
- double norm = ((riseTime + fallTime) / 2) * Math.sqrt(2 * Math.PI); //to ensure the total integral is equal to 1: = 33.8
-
- if (time < 3 * riseTime) {
- return corrGain * readoutGain * funcGaus(time - 3 * riseTime, riseTime) / norm;
- } else {
- return corrGain * readoutGain * funcGaus(time - 3 * riseTime, fallTime) / norm;
- }
+ double norm = ((riseTime + fallTime) / 2) * Math.sqrt(2 * Math.PI); //to ensure the total integral is equal to 1: = 33.8
+ return funcGaus(time - 3 * riseTime, (time < 3 * riseTime) ? riseTime : fallTime) / norm;
+ case ThreePole:
+ //peak at 2*tp
+ //peak value 2/(tp*e^2)
+ return ((time * time / (2 * shapingTime * shapingTime * shapingTime)) * Math.exp(-time / shapingTime));
+ default:
+ return 0.0;
}
}
java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal
--- java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal/FADCTriggerDriver.java 2014-04-21 19:03:07 UTC (rev 505)
+++ java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal/FADCTriggerDriver.java 2014-04-21 19:04:13 UTC (rev 506)
@@ -29,17 +29,18 @@
int nTriggers;
int totalEvents;
- protected double beamEnergy = 2.2 * ECalUtils.GeV;
+ protected double beamEnergy = -1; //by default, get beam energy from detector name
private int minHitCount = 1;
- private double clusterEnergyHigh = 1.85 / 2.2;
- private double clusterEnergyLow = .1 / 2.2;
- private double energySumThreshold = 1.0;
- private double energyDifferenceThreshold = 1.5 / 2.2;
- private double maxCoplanarityAngle = 35; // degrees
+ private boolean useDefaultCuts = true;
+ private double clusterEnergyHigh = 2.2 * ECalUtils.GeV;
+ private double clusterEnergyLow = .1 * ECalUtils.GeV;
+ private double energySumThreshold = 2.2 * ECalUtils.GeV;
+ private double energyDifferenceThreshold = 2.2 * ECalUtils.GeV;
+ private double maxCoplanarityAngle = 90; // degrees
// private double energyDistanceDistance = 250; // mm
// private double energyDistanceThreshold = 0.8 / 2.2;
private double energyDistanceDistance = 200; // mm
- private double energyDistanceThreshold = 0.5;
+ private double energyDistanceThreshold = 0.5; // unitless fraction
// maximum time difference between two clusters, in units of readout cycles (4 ns).
private int pairCoincidence = 2;
private double originX = 1393.0 * Math.tan(0.03052); //ECal midplane, defined by photon beam position (30.52 mrad) at ECal face (z=1393 mm)
@@ -84,27 +85,29 @@
this.clusterCollectionName = clusterCollectionName;
}
- public void setBeamEnergy(double beamEnergy) {
+ public void setCutsFromBeamEnergy(double beamEnergy) {
if (beamEnergy == 1.1) {
System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 1.1 GeV beam");
maxCoplanarityAngle = 90;
- clusterEnergyHigh = .7 / beamEnergy;
- clusterEnergyLow = .1 / beamEnergy;
- energySumThreshold = 0.8 / beamEnergy;
+ clusterEnergyHigh = .7 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 0.8 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
} else if (beamEnergy == 2.2) {
System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 2.2 GeV beam");
- maxCoplanarityAngle = 45;
- clusterEnergyHigh = 1.6 / beamEnergy;
- clusterEnergyLow = .1 / beamEnergy;
- energySumThreshold = 1.7 / beamEnergy;
+ maxCoplanarityAngle = 35;
+ clusterEnergyHigh = 1.5 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 1.9 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
} else if (beamEnergy == 6.6) {
System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 6.6 GeV beam");
maxCoplanarityAngle = 60;
- clusterEnergyHigh = 5.0 / beamEnergy;
- clusterEnergyLow = .1 / beamEnergy;
- energySumThreshold = 5.5 / beamEnergy;
+ clusterEnergyHigh = 5.0 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 5.5 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
}
- this.beamEnergy = beamEnergy * ECalUtils.GeV;
}
protected double getBeamEnergyFromDetector(Detector detector) {
@@ -129,7 +132,7 @@
/**
* Set X coordinate used as the origin for cluster coplanarity and distance
- * calculations. Defaults to the ECal midplane.
+ * calculations. Defaults to the ECal midplane. Units of mm.
*
* @param originX
*/
@@ -137,9 +140,107 @@
this.originX = originX;
}
+ /**
+ * Used for plot ranges and cuts that scale with energy. 1.1, 2.2 and 6.6
+ * are associated with default cuts. Units of GeV.
+ *
+ * @param beamEnergy
+ */
+ public void setBeamEnergy(double beamEnergy) {
+ this.beamEnergy = beamEnergy;
+ }
+
+ /**
+ * Use default cuts based on beam energy.
+ *
+ * @param useDefaultCuts
+ */
+ public void setUseDefaultCuts(boolean useDefaultCuts) {
+ this.useDefaultCuts = useDefaultCuts;
+ }
+
+ /**
+ * Minimum hit count for a cluster.
+ *
+ * @param minHitCount
+ */
+ public void setMinHitCount(int minHitCount) {
+ this.minHitCount = minHitCount;
+ }
+
+ /**
+ * Maximum energy for a cluster. Units of GeV.
+ *
+ * @param clusterEnergyHigh
+ */
+ public void setClusterEnergyHigh(double clusterEnergyHigh) {
+ this.clusterEnergyHigh = clusterEnergyHigh * ECalUtils.GeV;
+ }
+
+ /**
+ * Minimum energy for a cluster. Units of GeV.
+ *
+ * @param clusterEnergyLow
+ */
+ public void setClusterEnergyLow(double clusterEnergyLow) {
+ this.clusterEnergyLow = clusterEnergyLow * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum energy sum of the two clusters in a pair. Units of GeV.
+ *
+ * @param energySumThreshold
+ */
+ public void setEnergySumThreshold(double energySumThreshold) {
+ this.energySumThreshold = energySumThreshold * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum energy difference between the two clusters in a pair. Units of
+ * GeV.
+ *
+ * @param energyDifferenceThreshold
+ */
+ public void setEnergyDifferenceThreshold(double energyDifferenceThreshold) {
+ this.energyDifferenceThreshold = energyDifferenceThreshold * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum deviation from coplanarity for the two clusters in a pair. Units
+ * of degrees.
+ *
+ * @param maxCoplanarityAngle
+ */
+ public void setMaxCoplanarityAngle(double maxCoplanarityAngle) {
+ this.maxCoplanarityAngle = maxCoplanarityAngle;
+ }
+
+ /**
+ * Distance threshold for the energy-distance cut. Units of mm.
+ *
+ * @param energyDistanceDistance
+ */
+ public void setEnergyDistanceDistance(double energyDistanceDistance) {
+ this.energyDistanceDistance = energyDistanceDistance;
+ }
+
+ /**
+ * Energy threshold for the energy-distance cut. Units of the beam energy.
+ *
+ * @param energyDistanceThreshold
+ */
+ public void setEnergyDistanceThreshold(double energyDistanceThreshold) {
+ this.energyDistanceThreshold = energyDistanceThreshold;
+ }
+
@Override
public void detectorChanged(Detector detector) {
- setBeamEnergy(this.getBeamEnergyFromDetector(detector));
+ if (beamEnergy < 0) {
+ beamEnergy = this.getBeamEnergyFromDetector(detector);
+ }
+ if (useDefaultCuts) {
+ setCutsFromBeamEnergy(beamEnergy);
+ }
clusterHitCount2DAll = aida.histogram2D("All cluster pairs: hit count (less energetic vs. more energetic)", 9, 0.5, 9.5, 9, 0.5, 9.5);
clusterSumDiff2DAll = aida.histogram2D("All cluster pairs: energy difference vs. sum", 100, 0.0, 2 * beamEnergy, 100, 0.0, beamEnergy);
@@ -455,10 +556,10 @@
* @return true if a pair is found, false otherwise
*/
protected boolean clusterECut(HPSEcalCluster[] clusterPair) {
- return (clusterPair[0].getEnergy() < beamEnergy * clusterEnergyHigh
- && clusterPair[1].getEnergy() < beamEnergy * clusterEnergyHigh
- && clusterPair[0].getEnergy() > beamEnergy * clusterEnergyLow
- && clusterPair[1].getEnergy() > beamEnergy * clusterEnergyLow);
+ return (clusterPair[0].getEnergy() < clusterEnergyHigh
+ && clusterPair[1].getEnergy() < clusterEnergyHigh
+ && clusterPair[0].getEnergy() > clusterEnergyLow
+ && clusterPair[1].getEnergy() > clusterEnergyLow);
}
/**
@@ -470,7 +571,7 @@
*/
protected boolean energySum(Cluster[] clusterPair) {
double clusterESum = clusterPair[0].getEnergy() + clusterPair[1].getEnergy();
- return (clusterESum < beamEnergy * energySumThreshold);
+ return (clusterESum < energySumThreshold);
}
/**
@@ -483,7 +584,7 @@
protected boolean energyDifference(HPSEcalCluster[] clusterPair) {
double clusterEDifference = clusterPair[0].getEnergy() - clusterPair[1].getEnergy();
- return (clusterEDifference < beamEnergy * energyDifferenceThreshold);
+ return (clusterEDifference < energyDifferenceThreshold);
}
/**