Author: [log in to unmask]
Date: Wed Dec 7 13:48:50 2016
New Revision: 4609
Log:
final branch changes
Modified:
java/branches/branch-1116Fixes/ecal-recon/ (props changed)
java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverter.java
java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverterDriver.java
java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java
java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPOnlineClusterer.java
Modified: java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverter.java
=============================================================================
--- java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverter.java (original)
+++ java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverter.java Wed Dec 7 13:48:50 2016
@@ -1,328 +1,351 @@
-package org.hps.recon.ecal;
-
-import java.awt.event.ActionEvent;
-import java.awt.event.ActionListener;
-import java.util.ArrayList;
-import java.util.List;
-
-import org.hps.record.daqconfig.ConfigurationManager;
-import org.hps.record.daqconfig.FADCConfig;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.GenericObject;
-import org.lcsim.event.RawCalorimeterHit;
-import org.lcsim.event.RawTrackerHit;
-
-/**
- * <code>EcalOnlineRawConverter</code> handles the conversion of raw
- * hits of all modes to energy hit <code>CalorimeterHit</code> objects.
- * This converter will employ the runtime values for all parameters and
- * is intended to emulate the firmware specifically.<br/>
- * <br/>
- * The converter requires the presence of the DAQ configuration manager,
- * which is activated by either <code>DatabaseDAQConfigDriver</code>
- * or <code>DAQConfigDriver</code> depending on from where it is to
- * obtain the configuration.<br/>
- * <br/>
- * This converter is primarily employed in the trigger and hardware
- * diagnostic processes as well as the readout simulation in Monte
- * Carlo.
- *
- * @author Nathan Baltzell <[log in to unmask]>
- * @author Kyle McCarty <[log in to unmask]>
- */
-public class EcalOnlineRawConverter {
- // Defines the maximum number of peaks that may be extracted from
- // a single waveform.
- private int nPeak = 3;
- // The DAQ configuration manager for FADC parameters.
- private FADCConfig config = null;
- // Whether or not a constant integration window should be assumed
- // for the purpose of pedestal calculations/
- private boolean constantWindow = false;
- // The number of nanoseconds in a clock-cycle (sample).
- private static final int nsPerSample = 4;
-
- /**
- * Instantiates the <code>EcalOnlineRawConverter</code> and connects
- * it to the <code>ConfigurationManager</code> to receive settings
- * from the DAQ configuration.
- */
- public EcalOnlineRawConverter() {
- // Track changes in the DAQ configuration.
- ConfigurationManager.addActionListener(new ActionListener() {
- @Override
- public void actionPerformed(ActionEvent e) {
- // Get the FADC configuration.
- config = ConfigurationManager.getInstance().getFADCConfig();
-
- // Get the number of peaks.
- if(config.getMode() == 1) { nPeak = Integer.MAX_VALUE; }
- else { nPeak = config.getMaxPulses(); }
- }
- });
- }
-
- /**
- * Gets the pedestal for a given crystal and threshold crossing.
- * @param cellID - The cell ID of the crystal.
- * @param windowSamples - The size of the readout window. A value
- * of <code>-1</code> indicates an infinite readout window.
- * @param thresholdCrossing - The time of the threshold crossing in
- * 4-nanosecond clock-cycles (samples).
- * @return Returns the pedestal for the crystal and threshold
- * crossing.
- */
- public double getPulsePedestal(long cellID, int windowSamples, int thresholdCrossing) {
- // Track the starting and ending samples over which integration
- // will occur. Only the intermediary samples need be considered
- // for pedestal calculation.
- int firstSample, lastSample;
-
- // For finite readout windows, calculate the pedestal based on
- // the size of the full readout window in the event that the
- // integration window is larger than the readout window.
- if(windowSamples > 0 && (config.getNSA() + config.getNSB()) / nsPerSample >= windowSamples) {
- firstSample = 0;
- lastSample = windowSamples - 1;
- }
-
- // Otherwise, the pedestal should be calculated based on the
- // integration window size.
- else {
- // Define the sample width as equivalent to the integration
- // window size.
- firstSample = thresholdCrossing - config.getNSB() / nsPerSample;
- lastSample = thresholdCrossing + config.getNSA() / nsPerSample - 1;
-
- // In the event of a finite readout window, ignore any
- // samples that fall outside the readout window. Since these
- // are clipped and will not be integrated, these pedestals
- // do not contribute.
- if(windowSamples > 0) {
- if(firstSample < 0) { firstSample = 0; }
- if(lastSample >= windowSamples) { lastSample = windowSamples - 1; }
- }
- }
-
- // Calculate and return the pedestal.
- return(lastSample - firstSample + 1) * config.getPedestal(cellID);
- }
-
- /**
- * Converts a mode-1 digitized waveform into standard energy hit.
- * @param hit - The "hit" object representing the digitized waveform
- * for a given crystal.
- * @return Returns a list of <code>CalorimeterHit</code> objects
- * parsed from the waveform.
- */
- public List<CalorimeterHit> HitDtoA(RawTrackerHit hit) {
- // Get the cell ID for the crystal as well as the digitized
- // waveform samples.
- final long cellID = hit.getCellID();
- final short[] waveform = hit.getADCValues();
-
- // If there are no samples, then there is nothing to integrate.
- if(waveform.length == 0) { return null; }
-
- // The pulse integration threshold is defined as the combination
- // of the pedestal and the threshold configuration parameter.
- final int absoluteThreshold = (int) (config.getPedestal(cellID) + config.getThreshold(cellID));
-
- // Store each instance of a threshold crossing in that can be
- // found within the digitized waveform.
- List<Integer> thresholdCrossings = new ArrayList<Integer>();
-
- // Check for the special case of the first sample exceeding
- // the integration threshold.
- if(waveform[0] > absoluteThreshold) {
- thresholdCrossings.add(0);
- }
-
- // Search the remaining samples for threshold crossings.
- thresholdLoop:
- for(int sample = 1; sample < waveform.length; ++sample) {
- if(waveform[sample] > absoluteThreshold && waveform[sample - 1] <= absoluteThreshold) {
- // Add the sample index to the list of threshold crossing.
- thresholdCrossings.add(sample);
-
- // No new threshold crossings can be registered within
- // this pulse. In the case of mode-1 data, the end of
- // the pulse is considered to be 8 samples past the
- // crossing, as the per the SSP. Otherwise, it is defined
- // by the integration window.
- if(config.getMode() == 1) { sample += 8; }
- else { sample += config.getNSA() / nsPerSample - 1; }
-
- // If there is a limit defined on the maximum number
- // of peaks that may be processed, terminate the search
- // after this number of peaks have been found.
- if(thresholdCrossings.size() >= nPeak) { break thresholdLoop; }
- }
- }
-
- // Use the previously located threshold crossing to generate
- // calorimeter hits.
- List<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
- for(int thresholdCrossing : thresholdCrossings) {
- // Perform the pulse integral.
- final double[] data = convertWaveformToPulse(waveform, thresholdCrossing);
- double time = data[0];
- double sum = data[1];
-
- // Perform pedestal subtraction.
- sum -= getPulsePedestal(cellID, waveform.length, thresholdCrossing);
-
- // Perform gain scaling.
- double energy = adcToEnergy(sum, cellID);
-
- // Create a new hit and add it to the list.
- newHits.add(CalorimeterHitUtilities.create(energy, time, cellID));
- }
-
- // Return the list of hits.
- return newHits;
- }
-
- /**
- * Converts a raw mode-3 hit to a proper calorimeter hit in units
- * of energy.
- * @param hit - The raw hit that is to be converted.
- * @param timeOffset - The time offset for the hit.
- * @return Returns a <code>CalorimeterHit</code> hit object that
- * represents the raw mode-3 hit with units of energy and a correct
- * time-stamp.
- */
- public CalorimeterHit HitDtoA(RawCalorimeterHit hit, double timeOffset) {
- // Verify the validity of the time-stamp.
- if(hit.getTimeStamp() % 64 != 0) {
- System.out.println("Unexpected time-stamp: " + hit.getTimeStamp());
- }
-
- // Get the pedestal. In the case of a constant integration window
- // (i.e. infinite readout window size, so no risk of clipping),
- // the window width should be given as -1. Otherwise, the real
- // readout window size should be used.
- long id = hit.getCellID();
- double time = hit.getTimeStamp() / 16.0;
- double pedestal = getPulsePedestal(id, constantWindow ? -1 : config.getWindowWidth(), (int) time / nsPerSample);
-
- // Calculate the total ADC value for the pulse and convert it
- // to energy.
- double adcSum = hit.getAmplitude() - pedestal;
- double rawEnergy = adcToEnergy(adcSum, id);
-
- // Create a calorimeter hit from the result and return it.
- return CalorimeterHitUtilities.create(rawEnergy, time + timeOffset, id);
- }
-
- /**
- * Converts a raw mode-7 hit to a proper calorimeter hit in units
- * of energy.
- * @param hit - The raw hit that is to be converted.
- * @param mode7Data - Additional mode-7 data object.
- * @param timeOffset - The time offset for the hit.
- * @return Returns a <code>CalorimeterHit</code> hit object that
- * represents the raw mode-7 hit with units of energy and a correct
- * time-stamp.
- */
- public CalorimeterHit HitDtoA(RawCalorimeterHit hit, GenericObject mode7Data, double timeOffset) {
- // Get the time and crystal cell ID for the hit. Note that the
- // time-stamps use the full 62.5 ps resolution.
- double time = hit.getTimeStamp() / 16.0;
- long id = hit.getCellID();
-
- // Get the pedestal. In the case of a constant integration window
- // (i.e. infinite readout window size, so no risk of clipping),
- // the window width should be given as -1. Otherwise, the real
- // readout window size should be used.
- double pedestal = getPulsePedestal(id, constantWindow ? -1 : config.getWindowWidth(), (int) time / nsPerSample);
-
- // Calculate the total ADC value for the pulse and convert it
- // to energy.
- double adcSum = hit.getAmplitude() - pedestal;
- double rawEnergy = adcToEnergy(adcSum, id);
-
- // Create a calorimeter hit from the result and return it.
- return CalorimeterHitUtilities.create(rawEnergy, time + timeOffset, id);
- }
-
- /**
- * Converts a value in ADC in a crystal to energy in units of GeV.
- * @param adcSum - The ADC value to convert.
- * @param cellID - The cell ID of the crystal containing the value.
- * @return Returns the ADC value as an energy in units of GeV.
- */
- private double adcToEnergy(double adcSum, long cellID) {
- return config.getGain(cellID) * adcSum * EcalUtils.MeV;
- }
-
- /**
- * Emulate the FADC250 firmware in conversion of Mode-1 waveform to a Mode-3/7 pulse,
- * given a time for threshold crossing.
- */
-
- /**
- * Converts a mode-1 digitized waveform to a mode-3/7 pulse for a
- * a given threshold crossing.
- * @param waveform - The digitized waveform. Each array value should
- * correspond to a sample of the waveform in units of ADC.
- * @param thresholdCrossing - The time of the threshold crossing
- * in samples.
- * @return Returns a <code>double</code> primitive of size 2. The
- * first value represents the time in nanoseconds of the pulser and
- * the second value the total integrated value of the pulse in ADC.
- */
- private double[] convertWaveformToPulse(short[] waveform, int thresholdCrossing) {
- // Store the integration range.
- int firstSample, lastSample;
-
- // If the integration range is larger than the number of samples,
- // then all the samples are used for pulse integration.
- if((config.getNSA() + config.getNSB()) / nsPerSample >= waveform.length) {
- firstSample = 0;
- lastSample = waveform.length - 1;
- }
-
- // Otherwise, the integration range covers a number of samples
- // before and after the threshold crossing as defined by the
- // run parameters.
- else {
- firstSample = thresholdCrossing - config.getNSB() / nsPerSample;
- lastSample = thresholdCrossing + config.getNSA() / nsPerSample - 1;
- }
-
- // Perform the pulse integral.
- double sumADC = 0;
- integrationLoop:
- for (int sample = firstSample; sample <= lastSample; sample++) {
- // If the current sample occurs before the readout window,
- // then it does not exist and can not be integrated. Skip it.
- if(sample < 0) { continue integrationLoop; }
-
- // Likewise, samples that occur after the readout window are
- // not retained and must be skipped.
- if(sample >= waveform.length) { break integrationLoop; }
-
- // Otherwise, add the sample to the pulse total.
- sumADC += waveform[sample];
- }
-
- // Calculate the pulse time with a 4 nanosecond resolution.
- double pulseTime = thresholdCrossing * nsPerSample;
-
- // Return both the pulse time and the total integrated pulse ADC.
- return new double [] { pulseTime, sumADC };
- }
-
- /**
- * Sets whether to use a constant integration window for the the
- * purpose of determining the correct pedestal. This should be used
- * in conjunction with Monte Carlo data during the readout cycle.
- * @param state - <code>true</code> ignores the size of the readout
- * window when calculating pedestals, and <code>false</code> accounts
- * for it in the case of pulse-clipping.
- */
- void setUseConstantWindow(boolean state) {
- constantWindow = state;
- }
-}
+package org.hps.recon.ecal;
+
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.hps.record.daqconfig.ConfigurationManager;
+import org.hps.record.daqconfig.FADCConfig;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.GenericObject;
+import org.lcsim.event.RawCalorimeterHit;
+import org.lcsim.event.RawTrackerHit;
+
+/**
+ * <code>EcalOnlineRawConverter</code> handles the conversion of raw
+ * hits of all modes to energy hit <code>CalorimeterHit</code> objects.
+ * This converter will employ the runtime values for all parameters and
+ * is intended to emulate the firmware specifically.<br/>
+ * <br/>
+ * The converter requires the presence of the DAQ configuration manager,
+ * which is activated by either <code>DatabaseDAQConfigDriver</code>
+ * or <code>DAQConfigDriver</code> depending on from where it is to
+ * obtain the configuration.<br/>
+ * <br/>
+ * This converter is primarily employed in the trigger and hardware
+ * diagnostic processes as well as the readout simulation in Monte
+ * Carlo.
+ *
+ * @author Nathan Baltzell <[log in to unmask]>
+ * @author Kyle McCarty <[log in to unmask]>
+ */
+public class EcalOnlineRawConverter {
+ // Defines the maximum number of peaks that may be extracted from
+ // a single waveform.
+ private int nPeak = 3;
+ // The DAQ configuration manager for FADC parameters.
+ private FADCConfig config = null;
+ // Whether or not a constant integration window should be assumed
+ // for the purpose of pedestal calculations/
+ private boolean constantWindow = false;
+ // The number of nanoseconds in a clock-cycle (sample).
+ private static final int nsPerSample = 4;
+
+ /**
+ * Instantiates the <code>EcalOnlineRawConverter</code> and connects
+ * it to the <code>ConfigurationManager</code> to receive settings
+ * from the DAQ configuration.
+ */
+ public EcalOnlineRawConverter() {
+ // Track changes in the DAQ configuration.
+ ConfigurationManager.addActionListener(new ActionListener() {
+ @Override
+ public void actionPerformed(ActionEvent e) {
+ // Get the FADC configuration.
+ config = ConfigurationManager.getInstance().getFADCConfig();
+
+ // Get the number of peaks.
+ if(config.getMode() == 1) { nPeak = Integer.MAX_VALUE; }
+ else { nPeak = config.getMaxPulses(); }
+ }
+ });
+ }
+
+ /**
+ * Gets the pedestal for a given crystal and threshold crossing.
+ * @param cellID - The cell ID of the crystal.
+ * @param windowSamples - The size of the readout window. A value
+ * of <code>-1</code> indicates an infinite readout window.
+ * @param thresholdCrossing - The time of the threshold crossing in
+ * 4-nanosecond clock-cycles (samples).
+ * @return Returns the pedestal for the crystal and threshold
+ * crossing.
+ */
+ public int getPulsePedestal(long cellID, int windowSamples, int thresholdCrossing) {
+ // Track the starting and ending samples over which integration
+ // will occur. Only the intermediary samples need be considered
+ // for pedestal calculation.
+ int firstSample, lastSample;
+
+ // For finite readout windows, calculate the pedestal based on
+ // the size of the full readout window in the event that the
+ // integration window is larger than the readout window.
+ if(windowSamples > 0 && (config.getNSA() + config.getNSB()) / nsPerSample >= windowSamples) {
+ firstSample = 0;
+ lastSample = windowSamples - 1;
+ }
+
+ // Otherwise, the pedestal should be calculated based on the
+ // integration window size.
+ else {
+ // Define the sample width as equivalent to the integration
+ // window size.
+ firstSample = thresholdCrossing - config.getNSB() / nsPerSample;
+ lastSample = thresholdCrossing + config.getNSA() / nsPerSample - 1;
+
+ // In the event of a finite readout window, ignore any
+ // samples that fall outside the readout window. Since these
+ // are clipped and will not be integrated, these pedestals
+ // do not contribute.
+ if(windowSamples > 0) {
+ if(firstSample < 0) { firstSample = 0; }
+ if(lastSample >= windowSamples) { lastSample = windowSamples - 1; }
+ }
+ }
+
+ // Calculate and return the pedestal. The extra 1 MeV added to
+ // the pedestal offsets the rounding error that is incurred when
+ // converting it to an integer.
+ return (int) ((lastSample - firstSample + 1) * (config.getPedestal(cellID) + 0.001));
+ }
+
+ /**
+ * Converts a mode-1 digitized waveform into standard energy hit.
+ * @param hit - The "hit" object representing the digitized waveform
+ * for a given crystal.
+ * @return Returns a list of <code>CalorimeterHit</code> objects
+ * parsed from the waveform.
+ */
+ public List<CalorimeterHit> HitDtoA(RawTrackerHit hit) {
+ // Get the cell ID for the crystal as well as the digitized
+ // waveform samples.
+ final long cellID = hit.getCellID();
+ final short[] waveform = hit.getADCValues();
+
+ // If there are no samples, then there is nothing to integrate.
+ if(waveform.length == 0) { return null; }
+
+ // The pulse integration threshold is defined as the combination
+ // of the pedestal and the threshold configuration parameter.
+ final int absoluteThreshold = (int) (config.getPedestal(cellID) + config.getThreshold(cellID));
+
+ // Store each instance of a threshold crossing in that can be
+ // found within the digitized waveform.
+ List<Integer> thresholdCrossings = new ArrayList<Integer>();
+
+ // Check for the special case of the first sample exceeding
+ // the integration threshold.
+ if(waveform[0] > absoluteThreshold) {
+ thresholdCrossings.add(0);
+ }
+
+ // Search the remaining samples for threshold crossings.
+ thresholdLoop:
+ for(int sample = 1; sample < waveform.length; sample++) {
+ if(waveform[sample] > absoluteThreshold && waveform[sample - 1] <= absoluteThreshold) {
+ // Add the sample index to the list of threshold crossing.
+ thresholdCrossings.add(sample);
+
+ // No new threshold crossings can be registered within
+ // this pulse. In the case of mode-1 data, the end of
+ // the pulse is considered to be 8 samples past the
+ // crossing, as the per the SSP. Otherwise, it is defined
+ // by the integration window.
+ if(config.getMode() == 1) { sample += 8; }
+ else { sample += config.getNSA() / nsPerSample - 1; }
+
+ // If there is a limit defined on the maximum number
+ // of peaks that may be processed, terminate the search
+ // after this number of peaks have been found.
+ if(thresholdCrossings.size() >= nPeak) { break thresholdLoop; }
+ }
+ }
+
+ // Use the previously located threshold crossing to generate
+ // calorimeter hits.
+ List<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
+ for(int thresholdCrossing : thresholdCrossings) {
+ // Obtain the pedestal for the pulse.
+ int pedestal = getPulsePedestal(cellID, waveform.length, thresholdCrossing);
+
+ // Perform the pulse integral.
+ final int[] data = convertWaveformToPulse(waveform, thresholdCrossing);
+ int time = data[0];
+ int sum = data[1];
+
+ // Perform pedestal subtraction.
+ sum -= pedestal;
+
+ // Perform gain scaling.
+ double energy = adcToEnergy(sum, cellID);
+
+ // Hits should not have less than zero energy.
+ if(energy < 0) {
+ Logger.getGlobal().log(Level.WARNING, "A hit was produced with " + energy + " GeV energy!");
+ }
+
+ // Create a new hit and add it to the list.
+ newHits.add(CalorimeterHitUtilities.create(energy, time, cellID));
+ }
+
+ // Return the list of hits.
+ return newHits;
+ }
+
+ /**
+ * Converts a raw mode-3 hit to a proper calorimeter hit in units
+ * of energy.
+ * @param hit - The raw hit that is to be converted.
+ * @param timeOffset - The time offset for the hit.
+ * @return Returns a <code>CalorimeterHit</code> hit object that
+ * represents the raw mode-3 hit with units of energy and a correct
+ * time-stamp.
+ */
+ public CalorimeterHit HitDtoA(RawCalorimeterHit hit, double timeOffset) {
+ // Verify the validity of the time-stamp.
+ if(hit.getTimeStamp() % 64 != 0) {
+ System.out.println("Unexpected time-stamp: " + hit.getTimeStamp());
+ }
+
+ // Get the pedestal. In the case of a constant integration window
+ // (i.e. infinite readout window size, so no risk of clipping),
+ // the window width should be given as -1. Otherwise, the real
+ // readout window size should be used.
+ long id = hit.getCellID();
+ double time = hit.getTimeStamp() / 16.0;
+ int pedestal = getPulsePedestal(id, constantWindow ? -1 : config.getWindowWidth(), (int) time / nsPerSample);
+
+ // Calculate the total ADC value for the pulse and convert it
+ // to energy.
+ int adcSum = hit.getAmplitude() - pedestal;
+ double rawEnergy = adcToEnergy(adcSum, id);
+
+ // Create a calorimeter hit from the result and return it.
+ return CalorimeterHitUtilities.create(rawEnergy, time + timeOffset, id);
+ }
+
+ /**
+ * Converts a raw mode-7 hit to a proper calorimeter hit in units
+ * of energy.
+ * @param hit - The raw hit that is to be converted.
+ * @param mode7Data - Additional mode-7 data object.
+ * @param timeOffset - The time offset for the hit.
+ * @return Returns a <code>CalorimeterHit</code> hit object that
+ * represents the raw mode-7 hit with units of energy and a correct
+ * time-stamp.
+ */
+ public CalorimeterHit HitDtoA(RawCalorimeterHit hit, GenericObject mode7Data, double timeOffset) {
+ // Get the time and crystal cell ID for the hit. Note that the
+ // time-stamps use the full 62.5 ps resolution.
+ double time = hit.getTimeStamp() / 16.0;
+ long id = hit.getCellID();
+
+ // Get the pedestal. In the case of a constant integration window
+ // (i.e. infinite readout window size, so no risk of clipping),
+ // the window width should be given as -1. Otherwise, the real
+ // readout window size should be used.
+ int pedestal = getPulsePedestal(id, constantWindow ? -1 : config.getWindowWidth(), (int) time / nsPerSample);
+
+ // Calculate the total ADC value for the pulse and convert it
+ // to energy.
+ int adcSum = hit.getAmplitude() - pedestal;
+ double rawEnergy = adcToEnergy(adcSum, id);
+
+ // Create a calorimeter hit from the result and return it.
+ return CalorimeterHitUtilities.create(rawEnergy, time + timeOffset, id);
+ }
+
+ /**
+ * Converts a value in ADC in a crystal to energy in units of GeV.
+ * @param adcSum - The ADC value to convert.
+ * @param cellID - The cell ID of the crystal containing the value.
+ * @return Returns the ADC value as an energy in units of GeV.
+ */
+ private double adcToEnergy(int adcSum, long cellID) {
+ // Define the gain. To mimic the hardware, this is done through
+ // manipulation of integer values only. We multiply by 256 to
+ // preserve extra digits of accuracy.
+ int gain = (int) (256.0 * (config.getGain(cellID) + 0.001));
+
+ // Multiply the gain by the pulse-subtracted energy sum. Also
+ // remove the extra factor of 256. This gives the energy in units
+ // of MeV.
+ int energy = (int) ((adcSum * gain) / 256.0);
+
+ // Return the final energy as a double in units of GeV.
+ return energy * EcalUtils.MeV;
+ }
+
+ /**
+ * Emulate the FADC250 firmware in conversion of Mode-1 waveform to a Mode-3/7 pulse,
+ * given a time for threshold crossing.
+ */
+
+ /**
+ * Converts a mode-1 digitized waveform to a mode-3/7 pulse for a
+ * a given threshold crossing.
+ * @param waveform - The digitized waveform. Each array value should
+ * correspond to a sample of the waveform in units of ADC.
+ * @param thresholdCrossing - The time of the threshold crossing
+ * in samples.
+ * @return Returns a <code>double</code> primitive of size 2. The
+ * first value represents the time in nanoseconds of the pulser and
+ * the second value the total integrated value of the pulse in ADC.
+ */
+ private int[] convertWaveformToPulse(short[] waveform, int thresholdCrossing) {
+ // Store the integration range.
+ int firstSample, lastSample;
+
+ // If the integration range is larger than the number of samples,
+ // then all the samples are used for pulse integration.
+ if((config.getNSA() + config.getNSB()) / nsPerSample >= waveform.length) {
+ firstSample = 0;
+ lastSample = waveform.length - 1;
+ }
+
+ // Otherwise, the integration range covers a number of samples
+ // before and after the threshold crossing as defined by the
+ // run parameters.
+ else {
+ firstSample = thresholdCrossing - config.getNSB() / nsPerSample;
+ lastSample = thresholdCrossing + config.getNSA() / nsPerSample - 1;
+ }
+
+ // Perform the pulse integral.
+ int sumADC = 0;
+ integrationLoop:
+ for (int sample = firstSample; sample <= lastSample; sample++) {
+ // If the current sample occurs before the readout window,
+ // then it does not exist and can not be integrated. Skip it.
+ if(sample < 0) { continue integrationLoop; }
+
+ // Likewise, samples that occur after the readout window are
+ // not retained and must be skipped.
+ if(sample >= waveform.length) { break integrationLoop; }
+
+ // Otherwise, add the sample to the pulse total.
+ sumADC += waveform[sample];
+ }
+
+ // Calculate the pulse time with a 4 nanosecond resolution.
+ int pulseTime = thresholdCrossing * nsPerSample;
+
+ // Return both the pulse time and the total integrated pulse ADC.
+ return new int[] { pulseTime, sumADC };
+ }
+
+ /**
+ * Sets whether to use a constant integration window for the the
+ * purpose of determining the correct pedestal. This should be used
+ * in conjunction with Monte Carlo data during the readout cycle.
+ * @param state - <code>true</code> ignores the size of the readout
+ * window when calculating pedestals, and <code>false</code> accounts
+ * for it in the case of pulse-clipping.
+ */
+ void setUseConstantWindow(boolean state) {
+ constantWindow = state;
+ }
+}
Modified: java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverterDriver.java
=============================================================================
--- java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverterDriver.java (original)
+++ java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/EcalOnlineRawConverterDriver.java Wed Dec 7 13:48:50 2016
@@ -1,160 +1,160 @@
-package org.hps.recon.ecal;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.hps.record.daqconfig.ConfigurationManager;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.GenericObject;
-import org.lcsim.event.LCRelation;
-import org.lcsim.event.RawCalorimeterHit;
-import org.lcsim.event.RawTrackerHit;
-import org.lcsim.lcio.LCIOConstants;
-import org.lcsim.util.Driver;
-
-/**
- * This class is used to convert between collections of {@link org.lcsim.event.RawCalorimeterHit}
- * and {@link org.lcsim.event.RawTrackerHit}, objects with ADC/sample information, and
- * collections of {@link org.lcsim.event.CalorimeterHit}, objects with energy/time information.
- *
- * org.hps.recon.ecal.EcalRawConverter is called to do most of the lower level work.
- *
- *
-*/
-public class EcalOnlineRawConverterDriver extends Driver {
- private EcalOnlineRawConverter converter = null;
- /**
- * The input LCIO collection name. This can be either a
- * {@link org.lcsim.event.RawTrackerHit} or
- * {@link org.lcsim.event.RawCalorimeterHit}. These have ADC and
- * sample time information.
- */
- private String rawCollectionName = "EcalReadoutHits";
-
- /**
- * The output LCIO collection name. This will always a
- * {@link org.lcsim.event.CalorimeterHit} with energy (GeV) and
- * ns time information.
- */
- private String ecalCollectionName = "EcalCalHits";
-
- /**
- * ecalCollectionName "type" (must match detector-data)
- */
- private final String ecalReadoutName = "EcalHits";
-
- /**
- * Output relation between ecalCollectionName and Mode-7 pedestals.
- */
- private static final String extraDataRelationsName = "EcalReadoutExtraDataRelations";
-
- /**
- * Instantiates the <code>EcalOnlineRawConverter</code> for this
- * driver.
- */
- public EcalOnlineRawConverterDriver() { converter = new EcalOnlineRawConverter(); }
-
- /**
- * Checks that the required LCIO collection names are defined.
- */
- @Override
- public void startOfData() {
- if(ecalCollectionName == null) {
- throw new RuntimeException("The parameter ecalCollectionName was not set!");
- }
- }
-
- @Override
- public void process(EventHeader event) {
- // Do not process the event if the DAQ configuration is not
- // initialized. All online raw converter parameters are obtained
- // from this class, and this nothing can be done before they
- // are available.
- if(!ConfigurationManager.isInitialized()) {
- return;
- }
-
- double timeOffset = 0.0;
-
- // Define the LCIO data flags.
- int flags = 0;
- flags += 1 << LCIOConstants.RCHBIT_TIME; // Store hit time.
- flags += 1 << LCIOConstants.RCHBIT_LONG; // Store hit position; this flag has no effect for RawCalorimeterHits.
-
- // Create a list to store hits.
- ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
-
- // Events that contain RawTrackerHit objects use mode-1 data.
- if(event.hasCollection(RawTrackerHit.class, rawCollectionName)) {
- // Get the list of mode-1 waveforms.
- List<RawTrackerHit> hits = event.get(RawTrackerHit.class, rawCollectionName);
-
- // Extract hits from each waveform and store them.
- for(RawTrackerHit hit : hits) {
- newHits.addAll(converter.HitDtoA(hit));
- }
-
- // Add the hits to the data stream.
- event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName);
- }
-
- // Events that contain RawCalorimeterHit objects are either
- // mode-3 or mode-7.
- if(event.hasCollection(RawCalorimeterHit.class, rawCollectionName)) {
- // Check for extra relations data. This indicates that the
- // hits should be interpreted as mode-7.
- if(event.hasCollection(LCRelation.class, extraDataRelationsName)) { // extra information available from mode 7 readout
- List<LCRelation> extraDataRelations = event.get(LCRelation.class, extraDataRelationsName);
- for(LCRelation rel : extraDataRelations) {
- RawCalorimeterHit hit = (RawCalorimeterHit) rel.getFrom();
- GenericObject extraData = (GenericObject) rel.getTo();
- newHits.add(converter.HitDtoA(hit, extraData, timeOffset));
- }
- }
-
- // Otherwise, the hits should be treated as mode-3.
- else {
- List<RawCalorimeterHit> hits = event.get(RawCalorimeterHit.class, rawCollectionName);
- for(RawCalorimeterHit hit : hits) {
- newHits.add(converter.HitDtoA(hit, timeOffset));
- }
- }
- event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName);
- }
- }
-
- /**
- * Sets the output {@link org.lcsim.event.CalorimeterHit} LCIO
- * collection name.
- * @param ecalCollectionName - The LCIO collection name for output
- * data.
- */
- public void setEcalCollectionName(String ecalCollectionName) {
- this.ecalCollectionName = ecalCollectionName;
- }
-
- /**
- * Sets whether to use a constant integration window for the the
- * purpose of determining the correct pedestal. This should be used
- * in conjunction with Monte Carlo data during the readout cycle.
- * @param state - <code>true</code> ignores the size of the readout
- * window when calculating pedestals, and <code>false</code> accounts
- * for it in the case of pulse-clipping.
- */
- public void setIsReadoutSimulation(boolean state) {
- converter.setUseConstantWindow(state);
- }
-
- /**
- * Sets the input raw hit data LCIO collection name. Depending on
- * the driver configuration, this could be either a collection of
- * {@link org.lcsim.event.RawTrackerHit} objects for mode-1 data or
- * {@link org.lcsim.event.RawCalorimeterHit} objects for mode-3
- * and mode-7 data.
- * @param rawCollectionName - The LCIO collection name for raw data.
- */
- public void setRawCollectionName(String rawCollectionName) {
- this.rawCollectionName = rawCollectionName;
- }
-}
+package org.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.hps.record.daqconfig.ConfigurationManager;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.RawCalorimeterHit;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.lcio.LCIOConstants;
+import org.lcsim.util.Driver;
+
+/**
+ * This class is used to convert between collections of {@link org.lcsim.event.RawCalorimeterHit}
+ * and {@link org.lcsim.event.RawTrackerHit}, objects with ADC/sample information, and
+ * collections of {@link org.lcsim.event.CalorimeterHit}, objects with energy/time information.
+ *
+ * org.hps.recon.ecal.EcalRawConverter is called to do most of the lower level work.
+ *
+ *
+*/
+public class EcalOnlineRawConverterDriver extends Driver {
+ private EcalOnlineRawConverter converter = null;
+ /**
+ * The input LCIO collection name. This can be either a
+ * {@link org.lcsim.event.RawTrackerHit} or
+ * {@link org.lcsim.event.RawCalorimeterHit}. These have ADC and
+ * sample time information.
+ */
+ private String rawCollectionName = "EcalReadoutHits";
+
+ /**
+ * The output LCIO collection name. This will always a
+ * {@link org.lcsim.event.CalorimeterHit} with energy (GeV) and
+ * ns time information.
+ */
+ private String ecalCollectionName = "EcalCalHits";
+
+ /**
+ * ecalCollectionName "type" (must match detector-data)
+ */
+ private final String ecalReadoutName = "EcalHits";
+
+ /**
+ * Output relation between ecalCollectionName and Mode-7 pedestals.
+ */
+ private static final String extraDataRelationsName = "EcalReadoutExtraDataRelations";
+
+ /**
+ * Instantiates the <code>EcalOnlineRawConverter2</code> for this
+ * driver.
+ */
+ public EcalOnlineRawConverterDriver() { converter = new EcalOnlineRawConverter(); }
+
+ /**
+ * Checks that the required LCIO collection names are defined.
+ */
+ @Override
+ public void startOfData() {
+ if(ecalCollectionName == null) {
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+ }
+ }
+
+ @Override
+ public void process(EventHeader event) {
+ // Do not process the event if the DAQ configuration is not
+ // initialized. All online raw converter parameters are obtained
+ // from this class, and this nothing can be done before they
+ // are available.
+ if(!ConfigurationManager.isInitialized()) {
+ return;
+ }
+
+ double timeOffset = 0.0;
+
+ // Define the LCIO data flags.
+ int flags = 0;
+ flags += 1 << LCIOConstants.RCHBIT_TIME; // Store hit time.
+ flags += 1 << LCIOConstants.RCHBIT_LONG; // Store hit position; this flag has no effect for RawCalorimeterHits.
+
+ // Create a list to store hits.
+ ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
+
+ // Events that contain RawTrackerHit objects use mode-1 data.
+ if(event.hasCollection(RawTrackerHit.class, rawCollectionName)) {
+ // Get the list of mode-1 waveforms.
+ List<RawTrackerHit> hits = event.get(RawTrackerHit.class, rawCollectionName);
+
+ // Extract hits from each waveform and store them.
+ for(RawTrackerHit hit : hits) {
+ newHits.addAll(converter.HitDtoA(hit));
+ }
+
+ // Add the hits to the data stream.
+ event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName);
+ }
+
+ // Events that contain RawCalorimeterHit objects are either
+ // mode-3 or mode-7.
+ if(event.hasCollection(RawCalorimeterHit.class, rawCollectionName)) {
+ // Check for extra relations data. This indicates that the
+ // hits should be interpreted as mode-7.
+ if(event.hasCollection(LCRelation.class, extraDataRelationsName)) { // extra information available from mode 7 readout
+ List<LCRelation> extraDataRelations = event.get(LCRelation.class, extraDataRelationsName);
+ for(LCRelation rel : extraDataRelations) {
+ RawCalorimeterHit hit = (RawCalorimeterHit) rel.getFrom();
+ GenericObject extraData = (GenericObject) rel.getTo();
+ newHits.add(converter.HitDtoA(hit, extraData, timeOffset));
+ }
+ }
+
+ // Otherwise, the hits should be treated as mode-3.
+ else {
+ List<RawCalorimeterHit> hits = event.get(RawCalorimeterHit.class, rawCollectionName);
+ for(RawCalorimeterHit hit : hits) {
+ newHits.add(converter.HitDtoA(hit, timeOffset));
+ }
+ }
+ event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName);
+ }
+ }
+
+ /**
+ * Sets the output {@link org.lcsim.event.CalorimeterHit} LCIO
+ * collection name.
+ * @param ecalCollectionName - The LCIO collection name for output
+ * data.
+ */
+ public void setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ /**
+ * Sets whether to use a constant integration window for the the
+ * purpose of determining the correct pedestal. This should be used
+ * in conjunction with Monte Carlo data during the readout cycle.
+ * @param state - <code>true</code> ignores the size of the readout
+ * window when calculating pedestals, and <code>false</code> accounts
+ * for it in the case of pulse-clipping.
+ */
+ public void setIsReadoutSimulation(boolean state) {
+ converter.setUseConstantWindow(state);
+ }
+
+ /**
+ * Sets the input raw hit data LCIO collection name. Depending on
+ * the driver configuration, this could be either a collection of
+ * {@link org.lcsim.event.RawTrackerHit} objects for mode-1 data or
+ * {@link org.lcsim.event.RawCalorimeterHit} objects for mode-3
+ * and mode-7 data.
+ * @param rawCollectionName - The LCIO collection name for raw data.
+ */
+ public void setRawCollectionName(String rawCollectionName) {
+ this.rawCollectionName = rawCollectionName;
+ }
+}
Modified: java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java
=============================================================================
--- java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java (original)
+++ java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterEnergyCorrection.java Wed Dec 7 13:48:50 2016
@@ -1,6 +1,8 @@
package org.hps.recon.ecal.cluster;
import hep.physics.vec.Hep3Vector;
+
+import java.util.Random;
import org.hps.detector.ecal.EcalCrystal;
import org.hps.detector.ecal.HPSEcalDetectorElement;
@@ -19,27 +21,54 @@
* @author Jeremy McCormick <[log in to unmask]>
*/
public final class ClusterEnergyCorrection {
-
+
+
+ // Variables derived as the difference between data and mc noise in
+ // ecal cluster energy resolution.
+ static final double A = -0.0000981;
+ static final double B = 0.0013725;
+ static final double C = 0.00301;
+
+ // Calculate the noise factor to smear the Ecal energy by
+ private static double calcNoise(double energy){
+ Random r = new Random();
+ double noise = r.nextGaussian()*Math.sqrt(A+B*energy+C*Math.pow(energy, 2));
+ //System.out.println("energy:\t"+energy+"\tnoise:\t"+noise);
+ return noise;
+ }
+
+
// Variables for electron energy corrections.
static final double par0_em = -0.017;
- static final double par1_em[] = { 35, -0.06738, -0.0005613, 16.42, 0.3431,
- -2.021, 74.85, -0.3626 };
- static final double par2_em[] = { 35, 0.933, 0.003234, 18.06, 0.24, 8.586,
- 75.08, -0.39 };
+ static final double par1_em[] = { 35, -0.06738, -0.0005613, 16.42, 0.3431,-2.021, 74.85, -0.3626 };
+ static final double par2_em[] = { 35, 0.933, 0.003234, 18.06, 0.24, 8.586,75.08, -0.39 };
// Variables for positron energy corrections.
static final double par0_ep = -0.0131;
- static final double par1_ep[] = { 35, -0.076, -0.0008183, 17.88, 0.2886,
- -1.192, 73.12, -0.3747 };
- static final double par2_ep[] = { 35, 0.94, 0.003713, 18.19, 0.24, 8.342,
- 72.44, -0.39 };
+ static final double par1_ep[] = { 35, -0.076, -0.0008183, 17.88, 0.2886,-1.192, 73.12, -0.3747 };
+ static final double par2_ep[] = { 35, 0.94, 0.003713, 18.19, 0.24, 8.342,72.44, -0.39 };
// Variables for photon energy corrections.
static final double par0_p = -0.0113;
- static final double par1_p[] = { 35, -0.0585, -0.0008572, 16.76, 0.2784,
- -0.07232, 72.88, -0.1685 };
- static final double par2_p[] = { 35, 0.9307, 0.004, 18.05, 0.23, 3.027,
- 74.93, -0.34 };
+ static final double par1_p[] = { 35, -0.0585, -0.0008572, 16.76, 0.2784,-0.07232, 72.88, -0.1685 };
+ static final double par2_p[] = { 35, 0.9307, 0.004, 18.05, 0.23, 3.027,74.93, -0.34 };
+
+ // Variables for electron energy corrections--MC.
+ static final double par0MC_em = 0.009051;
+ static final double par1MC_em[] = {35,-0.1322,-0.0005613,16.42,0.3431,-2.021,74.85,-0.3626};
+ static final double par2MC_em[] = {35, 0.9652, 0.003234, 18.06, 0.2592, 8.586, 75.08, -0.3771};
+
+
+
+ // Variables for positron energy corrections--MC.
+ static final double par0MC_ep = 0.01307;
+ static final double par1MC_ep[] = {35,-0.1415,-0.0008183,17.88,0.2886,-1.192,73.12,-0.3747};
+ static final double par2MC_ep[] = {35, 0.9733, 0.003713, 18.19, 0.2557, 8.342, 72.44, -0.3834};
+
+ // Variables for photon energy corrections--MC.
+ static final double par0MC_p = 0.01604;
+ static final double par1MC_p[] = {35,-0.1268,-0.0008572,16.76,0.2784,-0.07232,72.88,-0.1685};
+ static final double par2MC_p[] = {35, 0.965, 0.004, 18.05, 0.24, 3.027, 74.93, -0.3221};
/**
* Calculate the corrected energy for the cluster.
@@ -48,10 +77,9 @@
* The input cluster.
* @return The corrected energy.
*/
- public static double calculateCorrectedEnergy(HPSEcal3 ecal, Cluster cluster) {
+ public static double calculateCorrectedEnergy(HPSEcal3 ecal, Cluster cluster,boolean isMC) {
double rawE = cluster.getEnergy();
- return computeCorrectedEnergy(ecal, cluster.getParticleId(), rawE,
- cluster.getPosition()[0], cluster.getPosition()[1]);
+ return computeCorrectedEnergy(ecal, cluster.getParticleId(), rawE,cluster.getPosition()[0], cluster.getPosition()[1],isMC);
}
/**
@@ -62,11 +90,9 @@
* The input cluster.
* @return The corrected energy.
*/
- public static double calculateCorrectedEnergy(HPSEcal3 ecal,
- Cluster cluster, double ypos) {
+ public static double calculateCorrectedEnergy(HPSEcal3 ecal,Cluster cluster, double ypos, boolean isMC) {
double rawE = cluster.getEnergy();
- return computeCorrectedEnergy(ecal, cluster.getParticleId(), rawE,
- cluster.getPosition()[0], ypos);
+ return computeCorrectedEnergy(ecal, cluster.getParticleId(), rawE, cluster.getPosition()[0], ypos,isMC);
}
/**
@@ -75,8 +101,11 @@
* @param cluster
* The input cluster.
*/
- public static void setCorrectedEnergy(HPSEcal3 ecal, BaseCluster cluster) {
- double correctedEnergy = calculateCorrectedEnergy(ecal, cluster);
+ public static void setCorrectedEnergy(HPSEcal3 ecal, BaseCluster cluster, boolean isMC) {
+ double correctedEnergy = calculateCorrectedEnergy(ecal, cluster,isMC);
+ if (isMC){
+ correctedEnergy += calcNoise(correctedEnergy);
+ }
cluster.setEnergy(correctedEnergy);
}
@@ -87,9 +116,11 @@
* The input cluster.
*/
- public static void setCorrectedEnergy(HPSEcal3 ecal, BaseCluster cluster,
- double ypos) {
- double correctedEnergy = calculateCorrectedEnergy(ecal, cluster, ypos);
+ public static void setCorrectedEnergy(HPSEcal3 ecal, BaseCluster cluster,double ypos, boolean isMC) {
+ double correctedEnergy = calculateCorrectedEnergy(ecal, cluster, ypos,isMC);
+ if(isMC){
+ correctedEnergy += calcNoise(correctedEnergy);
+ }
cluster.setEnergy(correctedEnergy);
}
@@ -107,8 +138,7 @@
* @return Corrected Energy
*/
- private static double computeCorrectedEnergy(HPSEcal3 ecal, int pdg,
- double rawEnergy, double xpos, double ypos) {
+ private static double computeCorrectedEnergy(HPSEcal3 ecal, int pdg, double rawEnergy, double xpos, double ypos, boolean isMC) {
// distance to beam gap edge
double r;
// Get these values from the Ecal geometry:
@@ -118,24 +148,20 @@
// 22.3;//ecal.getNode().getChild("layout").getAttribute("beamgapTop").getDoubleValue();//mm
double BEAMGAPTOP = 20.0;
try {
- BEAMGAPTOP = ecal.getNode().getChild("layout")
- .getAttribute("beamgapTop").getDoubleValue();
+ BEAMGAPTOP = ecal.getNode().getChild("layout").getAttribute("beamgapTop").getDoubleValue();
} catch (Exception e) {
try {
- BEAMGAPTOP = ecal.getNode().getChild("layout")
- .getAttribute("beamgap").getDoubleValue();
+ BEAMGAPTOP = ecal.getNode().getChild("layout").getAttribute("beamgap").getDoubleValue();
} catch (Exception ee) {
ee.printStackTrace();
}
}
double BEAMGAPBOT = -20.0;
try {
- BEAMGAPBOT = -ecal.getNode().getChild("layout")
- .getAttribute("beamgapBottom").getDoubleValue();
+ BEAMGAPBOT = -ecal.getNode().getChild("layout").getAttribute("beamgapBottom").getDoubleValue();
} catch (Exception e) {
try {
- BEAMGAPBOT = ecal.getNode().getChild("layout")
- .getAttribute("beamgap").getDoubleValue();
+ BEAMGAPBOT = ecal.getNode().getChild("layout").getAttribute("beamgap").getDoubleValue();
} catch (Exception ee) {
ee.printStackTrace();
}
@@ -175,24 +201,41 @@
//Eliminates corrections at outermost edges to negative cluster energies
//66 for positrons, 69 is safe for electrons and photons
if (r > 66) {r = 66;}
-
- switch (pdg) {
- case 11:
- // electron
- return computeCorrectedEnergy(r, rawEnergy, par0_em, par1_em,
- par2_em);
- case -11:
- // positron
- return computeCorrectedEnergy(r, rawEnergy, par0_ep, par1_ep,
- par2_ep);
- case 22:
- // photon
- return computeCorrectedEnergy(r, rawEnergy, par0_p, par1_p, par2_p);
- default:
- // unknown
- return rawEnergy;
- }
- }
+
+ if (isMC){
+ switch (pdg) {
+ case 11:
+ // electron
+ return computeCorrectedEnergy(r, rawEnergy, par0MC_em, par1MC_em,par2MC_em);
+ case -11:
+ // positron
+ return computeCorrectedEnergy(r, rawEnergy, par0MC_ep, par1MC_ep,par2MC_ep);
+ case 22:
+ // photon
+ return computeCorrectedEnergy(r, rawEnergy, par0MC_p, par1MC_p, par2MC_p);
+ default:
+ // unknown
+ return rawEnergy;
+ }
+ }
+ else{
+ switch (pdg) {
+ case 11:
+ // electron
+ return computeCorrectedEnergy(r, rawEnergy, par0_em, par1_em,par2_em);
+ case -11:
+ // positron
+ return computeCorrectedEnergy(r, rawEnergy, par0_ep, par1_ep,par2_ep);
+ case 22:
+ // photon
+ return computeCorrectedEnergy(r, rawEnergy, par0_p, par1_p, par2_p);
+ default:
+ // unknown
+ return rawEnergy;
+ }
+ }
+ }
+
/**
* Calculates the energy correction to a cluster given the variables from
Modified: java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPOnlineClusterer.java
=============================================================================
--- java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPOnlineClusterer.java (original)
+++ java/branches/branch-1116Fixes/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPOnlineClusterer.java Wed Dec 7 13:48:50 2016
@@ -108,7 +108,6 @@
private IHistogram1D clusterTotalEnergy = aida.histogram1D("GTP(O) Cluster Plots/Cluster Total Energy Distribution", 176, 0.0, 2.2);
private IHistogram2D hitDistribution = aida.histogram2D("GTP(O) Cluster Plots/Hit Distribution", 46, -23, 23, 11, -5.5, 5.5);
private IHistogram2D clusterDistribution = aida.histogram2D("GTP(O) Cluster Plots/Cluster Seed Distribution", 46, -23, 23, 11, -5.5, 5.5);
- private IHistogram1D energyDistribution = aida.histogram1D("GTP(O) Cluster Plots/Percent Negative Energy Distribution", 100, 0.0, 1.0);
/**
* Instantiates a new instance of a readout GTP clustering algorithm.
@@ -140,10 +139,10 @@
Collections.sort(hitList, new Comparator<CalorimeterHit>() {
@Override
public int compare(CalorimeterHit firstHit, CalorimeterHit secondHit) {
- int[] ix = { firstHit.getIdentifierFieldValue("ix"), secondHit.getIdentifierFieldValue("ix") };
+ int[] ix = { getHitX(firstHit), getHitX(secondHit) };
if(ix[0] != ix[1]) { return Integer.compare(ix[0], ix[1]); }
else {
- int iy[] = { firstHit.getIdentifierFieldValue("iy"), secondHit.getIdentifierFieldValue("iy") };
+ int iy[] = { getHitY(firstHit), getHitY(secondHit) };
return Integer.compare(iy[0], iy[1]);
}
}
@@ -152,12 +151,7 @@
// Print the hit collection.
System.out.println("Event Hit Collection:");
for(CalorimeterHit hit : hitList) {
- int ix = hit.getIdentifierFieldValue("ix");
- int iy = hit.getIdentifierFieldValue("iy");
- double energy = hit.getCorrectedEnergy();
- double time = hit.getTime();
-
- System.out.printf("\tHit --> %6.3f GeV at (%3d, %3d) and at t = %.2f%n", energy, ix, iy, time);
+ System.out.printf("\t%s%n", getHitText(hit));
}
System.out.println();
}
@@ -169,7 +163,7 @@
Collections.sort(hitList, new Comparator<CalorimeterHit>() {
@Override
public int compare(CalorimeterHit firstHit, CalorimeterHit secondHit) {
- return Double.compare(secondHit.getTime(), firstHit.getTime());
+ return Double.compare(getHitTime(secondHit), getHitTime(firstHit));
}
});
@@ -183,15 +177,24 @@
// Iterate over each hit and see if it qualifies as a seed hit.
seedLoop:
for(CalorimeterHit seed : hitList) {
+ // VERBOSE :: Output the seed that is being considered.
+ if(verbose) {
+ System.out.println("\n");
+ System.out.println("Considering seed " + getHitText(seed));
+ }
+
// Put the hit energy into the hit energy distribution.
- hitEnergy.fill(seed.getCorrectedEnergy());
- hitDistribution.fill(seed.getIdentifierFieldValue("ix"), seed.getIdentifierFieldValue("iy"));
+ hitEnergy.fill(getHitEnergy(seed));
+ hitDistribution.fill(getHitX(seed), getHitY(seed));
// Check whether the potential seed passes the seed
// energy cut.
- if(seed.getCorrectedEnergy() < seedThreshold) {
+ if(verbose) { System.out.printf("Checking seed energy threshold %5.3f >= %5.3f... ", getHitEnergy(seed), seedThreshold); }
+ if(getHitEnergy(seed) < seedThreshold) {
+ if(verbose) { System.out.println("[fail]"); }
continue seedLoop;
}
+ if(verbose) { System.out.println("[pass]"); }
// Create a cluster for the potential seed.
BaseCluster protoCluster = createBasicCluster();
@@ -247,16 +250,6 @@
clusterHitCount.fill(protoCluster.getCalorimeterHits().size());
clusterDistribution.fill(protoCluster.getCalorimeterHits().get(0).getIdentifierFieldValue("ix"),
protoCluster.getCalorimeterHits().get(0).getIdentifierFieldValue("iy"));
-
- // Determine how much energy in the cluster is negative
- // and how is positive.
- double nenergy = 0.0;
- double penergy = 0.0;
- for(CalorimeterHit hit : protoCluster.getCalorimeterHits()) {
- if(hit.getCorrectedEnergy() > 0) { penergy += hit.getCorrectedEnergy(); }
- else { nenergy += hit.getCorrectedEnergy(); }
- }
- energyDistribution.fill(Math.abs(nenergy) / (penergy + Math.abs(nenergy)));
}
// VERBOSE :: Print out all the clusters in the event.
@@ -391,6 +384,26 @@
// windows appropriately.
timeAfter = cyclesAfter * 4;
timeWindow = Math.max(timeBefore, timeAfter);
+ }
+
+ private static final String getHitText(CalorimeterHit hit) {
+ return String.format("Hit --> %6.3f GeV at (%3d, %3d) and at t = %.2f", getHitEnergy(hit), getHitX(hit), getHitY(hit), getHitTime(hit));
+ }
+
+ private static final double getHitEnergy(CalorimeterHit hit) {
+ return hit.getCorrectedEnergy();
+ }
+
+ private static final int getHitX(CalorimeterHit hit) {
+ return hit.getIdentifierFieldValue("ix");
+ }
+
+ private static final int getHitY(CalorimeterHit hit) {
+ return hit.getIdentifierFieldValue("iy");
+ }
+
+ private static final double getHitTime(CalorimeterHit hit) {
+ return hit.getTime();
}
/**
@@ -561,4 +574,4 @@
// treated as within time.
else { return false; }
}
-}
+}
|