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-24 12:08:22 UTC (rev 509)
+++ java/branches/ecal-readout-sim_HPSJAVA-93/src/main/java/org/hps/readout/ecal/FADCEcalReadoutDriver.java 2014-04-24 13:04:56 UTC (rev 510)
@@ -15,7 +15,12 @@
import java.util.PriorityQueue;
import java.util.Set;
-import org.hps.conditions.deprecated.EcalConditions;
+import org.hps.conditions.TableConstants;
+import org.hps.conditions.ecal.EcalChannel.EcalChannelCollection;
+import org.hps.conditions.ecal.EcalChannel.GeometryId;
+import org.hps.conditions.ecal.EcalChannelConstants;
+import org.hps.conditions.ecal.EcalConditions;
+import org.lcsim.conditions.ConditionsManager;
import org.hps.recon.ecal.ECalUtils;
import org.hps.recon.ecal.HPSRawCalorimeterHit;
import org.lcsim.event.CalorimeterHit;
@@ -29,7 +34,13 @@
import org.lcsim.geometry.subdetector.HPSEcal3;
import org.hps.util.RandomGaussian;
import org.lcsim.lcio.LCIOConstants;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.detector.identifier.Identifier;
+
+
+
/**
* Performs readout of ECal hits. Simulates time evolution of preamp output
* pulse.
@@ -45,6 +56,11 @@
private static final int ECAL_PULSE_INTEGRAL_MODE = 3;
String ecalName = "Ecal";
Subdetector ecal;
+
+ EcalConditions ecalConditions = null;
+ IIdentifierHelper helper = null;
+ EcalChannelCollection channels = null;
+
//buffer for preamp signals (units of volts, no pedestal)
private Map<Long, RingBuffer> signalMap = null;
//ADC pipeline for readout (units of ADC counts)
@@ -83,9 +99,14 @@
//output collection name for hits read out from trigger
private String ecalReadoutCollectionName = "EcalReadoutHits";
private int mode = ECAL_PULSE_INTEGRAL_MODE;
+
private int readoutThreshold = 10;
private int triggerThreshold = 10;
- private double scaleFactor = 1;
+ //amplitude ADC counts/GeV
+// private double gain = 0.5*1000 * 80.0 / 60;
+ private double scaleFactor = 128;
+
+
private double fixedGain = -1;
private boolean constantTriggerWindow = true;
private boolean addNoise = false;
@@ -220,21 +241,33 @@
public Map<Long, FADCPipeline> getPipelineMap() {
return pipelineMap;
}
-
@Override
protected void readHits(List<RawCalorimeterHit> hits) {
for (Long cellID : signalMap.keySet()) {
+
RingBuffer signalBuffer = signalMap.get(cellID);
-
FADCPipeline pipeline = pipelineMap.get(cellID);
pipeline.step();
+
+ // Get the channel data.
+ EcalChannelConstants channelData = findChannel(cellID);
double currentValue = signalBuffer.currentValue() * ((Math.pow(2, nBit) - 1) / maxVolt); //12-bit ADC with maxVolt V range
- int pedestal = (int) Math.round(EcalConditions.physicalToPedestal(cellID));
+
+ /*
+ double pedestal = channelData.getCalibration().getPedestal();
+ 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
+*/
+
+// Pedestals will likely be expressed in number of fADC counts, so they should be int.
+ int pedestal = (int) channelData.getCalibration().getPedestal();
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
+
+// I think that it is easier to read the code when digitizedValue is present, is there any coding reason not to do so?
pipeline.writeValue(digitizedValue);
int pedestalSubtractedValue = digitizedValue - pedestal;
+
//System.out.println(signalBuffer.currentValue() + " " + currentValue + " " + pipeline.currentValue());
Integer sum = sumMap.get(cellID);
@@ -369,6 +402,10 @@
// System.out.println("Reading FADC data");
List<RawTrackerHit> hits = new ArrayList<RawTrackerHit>();
for (Long cellID : pipelineMap.keySet()) {
+
+ // Get the channel data.
+ EcalChannelConstants channelData = findChannel(cellID);
+
short[] window = getWindow(cellID);
short[] adcValues = null;
int pointerOffset = 0;
@@ -381,7 +418,7 @@
if (numSamplesToRead == 0) {
hits.add(new BaseRawTrackerHit(cellID, thresholdCrossing, adcValues));
}
- } else if ((i == 0 || window[i - 1] <= EcalConditions.physicalToPedestal(cellID) + readoutThreshold) && window[i] > EcalConditions.physicalToPedestal(cellID) + readoutThreshold) {
+ } else if ((i == 0 || window[i - 1] <= channelData.getCalibration().getPedestal() + readoutThreshold) && window[i] > channelData.getCalibration().getPedestal() + readoutThreshold) {
thresholdCrossing = i;
pointerOffset = Math.min(numSamplesBefore, i);
numSamplesToRead = pointerOffset + Math.min(numSamplesAfter, readoutWindow - i - pointerOffset - 1);
@@ -396,6 +433,10 @@
// System.out.println("Reading FADC data");
List<RawCalorimeterHit> hits = new ArrayList<RawCalorimeterHit>();
for (Long cellID : pipelineMap.keySet()) {
+
+ // Get the channel data.
+ EcalChannelConstants channelData = findChannel(cellID);
+
short[] window = getWindow(cellID);
int adcSum = 0;
int pointerOffset = 0;
@@ -412,7 +453,7 @@
if (numSamplesToRead == 0) {
hits.add(new BaseRawCalorimeterHit(cellID, adcSum, 64 * thresholdCrossing));
}
- } else if ((i == 0 || window[i - 1] <= EcalConditions.physicalToPedestal(cellID) + readoutThreshold) && window[i] > EcalConditions.physicalToPedestal(cellID) + readoutThreshold) {
+ } else if ((i == 0 || window[i - 1] <= channelData.getCalibration().getPedestal() + readoutThreshold) && window[i] > channelData.getCalibration().getPedestal() + readoutThreshold) {
thresholdCrossing = i;
pointerOffset = Math.min(numSamplesBefore, i);
numSamplesToRead = pointerOffset + Math.min(numSamplesAfter, readoutWindow - i - pointerOffset - 1);
@@ -430,13 +471,21 @@
for (CalorimeterHit hit : hits) {
RingBuffer eDepBuffer = signalMap.get(hit.getCellID());
double energyAmplitude = hit.getRawEnergy();
+ // Get the channel data.
+ EcalChannelConstants channelData = findChannel(hit.getCellID());
+
if (addNoise) {
//add preamp noise and photoelectron Poisson noise in quadrature
double noise;
+<<<<<<< .mine
+ if (!useCRRCShape) {
+ noise = Math.sqrt(Math.pow(channelData.getCalibration().getNoise() * channelData.getGain().getGain() * 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));
+>>>>>>> .r508
} else {
- noise = Math.sqrt(Math.pow(EcalConditions.physicalToNoise(hit.getCellID()) * EcalConditions.physicalToGain(hit.getCellID()) * ECalUtils.MeV, 2) + hit.getRawEnergy() * ECalUtils.MeV / pePerMeV);
+ noise = Math.sqrt(Math.pow(channelData.getCalibration().getNoise() * channelData.getGain().getGain() * ECalUtils.MeV, 2) + hit.getRawEnergy() * ECalUtils.MeV / pePerMeV);
}
energyAmplitude += RandomGaussian.getGaussian(0, noise);
}
@@ -458,7 +507,18 @@
@Override
public void detectorChanged(Detector detector) {
// Get the Subdetector.
- ecal = detector.getSubdetector(ecalName);
+ //ecal = detector.getSubdetector(ecalName);
+
+ // ECAL combined conditions object.
+ ecalConditions = ConditionsManager.defaultInstance()
+ .getCachedConditions(EcalConditions.class, TableConstants.ECAL_CONDITIONS).getCachedData();
+
+ // List of channels.
+ channels = ecalConditions.getChannelMap();
+
+ // ID helper.
+ helper = detector.getSubdetector("Ecal").getDetectorElement().getIdentifierHelper();
+
resetFADCBuffers();
}
@@ -470,13 +530,23 @@
pipelineMap = new HashMap<Long, FADCPipeline>();
Set<Long> cells = ((HPSEcal3) ecal).getNeighborMap().keySet();
for (Long cellID : cells) {
+ EcalChannelConstants channelData = findChannel(cellID);
signalMap.put(cellID, new RingBuffer(bufferLength));
- pipelineMap.put(cellID, new FADCPipeline(pipelineLength, (int) Math.round(EcalConditions.physicalToPedestal(cellID))));
+ pipelineMap.put(cellID, new FADCPipeline(pipelineLength, (int) Math.round(channelData.getCalibration().getPedestal())));
}
return true;
}
private double pulseAmplitude(double time, long cellID) {
+<<<<<<< .mine
+
+ // Get the channel data.
+ EcalChannelConstants channelData = findChannel(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;
@@ -484,8 +554,8 @@
corrGain = fixedGain;
} else {
corrGain = 1.0 / EcalConditions.physicalToGain(cellID);
+>>>>>>> .r508
}
-
return corrGain * readoutGain * pulseAmplitude(time, pulseShape, tp);
} else {
//normalization constant from cal gain (MeV/integral bit) to amplitude gain (amplitude bit/GeV)
@@ -493,13 +563,22 @@
if (fixedGain > 0) {
gain = readoutPeriod / (fixedGain * ECalUtils.MeV * ((Math.pow(2, nBit) - 1) / maxVolt));
} else {
- gain = readoutPeriod / (EcalConditions.physicalToGain(cellID) * ECalUtils.MeV * ((Math.pow(2, nBit) - 1) / maxVolt));
+ gain = readoutPeriod / (channelData.getGain().getGain() * ECalUtils.MeV * ((Math.pow(2, nBit) - 1) / maxVolt));
}
return gain * pulseAmplitude(time, pulseShape, tp);
}
}
+<<<<<<< .mine
+ //if fixedGain is set, multiply the default gain by this factor
+ double corrGain = 1.0;
+ if (fixedGain > 0) {
+ corrGain = fixedGain;
+ } else {
+ corrGain = 1.0 / channelData.getGain().getGain();
+ }
+=======
/**
* Returns pulse amplitude at the given time (relative to hit time).
* Amplitude is normalized so the pulse integral is 1.
@@ -520,6 +599,7 @@
//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
+>>>>>>> .r508
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;
@@ -584,4 +664,27 @@
return array[((ptr - pos) % size + size) % size];
}
}
+
+ // Convert physical ID to gain value.
+ private EcalChannelConstants findChannel(long cellID) {
+ // Make an ID object from raw hit ID.
+ IIdentifier id = new Identifier(cellID);
+
+ // Get physical field values.
+ int x = helper.getValue(id, "ix");
+ int y = helper.getValue(id, "iy");
+
+ // Create an ID to search for in channel collection.
+ GeometryId geometryId = new GeometryId();
+ geometryId.x = x;
+ geometryId.y = y;
+
+ // Find the ECAL channel.
+// return channels.findChannel(geometryId);
+
+ // Get the channel data.
+ return ecalConditions.getChannelConstants(channels.findChannel(geometryId));
+ }
+
+
}
