diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
index 7b9816c36241ca1b9a4041c046f0319ed79070a8..7671d13755f37e945b679796a15cb69f5d7aba98 100755
--- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
+++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py
@@ -24,6 +24,8 @@ parser.add_argument("-r", "--reco", default="",
help="Specify reco tag (to append to output filename)")
parser.add_argument("-n", "--nevents", type=int, default=-1,
help="Specify number of events to process (default: all)")
+parser.add_argument("-v", "--verbose", action='store_true',
+ help="Turn on DEBUG output")
args = parser.parse_args()
@@ -165,17 +167,20 @@ replicaSvc.UseGeomSQLite = True
# Configure verbosity
# ConfigFlags.dump()
-# logging.getLogger('forcomps').setLevel(VERBOSE)
-# acc.foreach_component("*").OutputLevel = VERBOSE
-acc.foreach_component("*").OutputLevel = INFO
+if args.verbose:
+ acc.foreach_component("*").OutputLevel = VERBOSE
+
+ #acc.getService("FaserByteStreamInputSvc").DumpFlag = True
+ #acc.getService("FaserEventSelector").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE
+ #acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE
+
+else:
+ acc.foreach_component("*").OutputLevel = INFO
+
acc.foreach_component("*ClassID*").OutputLevel = INFO
-# log.setLevel(VERBOSE)
-#acc.getService("FaserByteStreamInputSvc").DumpFlag = True
-#acc.getService("FaserEventSelector").OutputLevel = VERBOSE
-#acc.getService("FaserByteStreamInputSvc").OutputLevel = VERBOSE
-#acc.getService("FaserByteStreamCnvSvc").OutputLevel = VERBOSE
-#acc.getService("FaserByteStreamAddressProviderSvc").OutputLevel = VERBOSE
acc.getService("MessageSvc").Format = "% F%40W%S%7W%R%T %0W%M"
# Execute and finish
diff --git a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
index 3c6da1c0eb5a365f469dba9acd0b2027c433c2f4..f1320c41e7dc75f5c9c8cd701dc06a92c3bab1c3 100644
--- a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
+++ b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
@@ -69,16 +69,8 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' initialized");
- // Reconstruct each waveform
- for( const auto& wave : *waveformHandle) {
-
- ATH_MSG_DEBUG("Reconstruct waveform for channel " << wave->channel());
-
- // Reconstruct the hits, may be more than one, so pass container
- CHECK( m_recoTool->reconstruct(*wave, clockptr,
- hitContainerHandle.ptr()) );
-
- }
+ // Reconstruct all waveforms
+ CHECK( m_recoTool->reconstructAll(*waveformHandle, clockptr, hitContainerHandle.ptr()) );
ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' filled with "<< hitContainerHandle->size() <<" items");
diff --git a/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h b/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
index e985cb1d66b81e0bddbee05e682798fdb0f61be8..cc10197b262f0695327f7304a84998d198300678 100644
--- a/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
+++ b/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
@@ -20,6 +20,7 @@
#include "xAODFaserWaveform/WaveformClock.h"
class RawWaveform;
+class RawWaveformContainer;
///Interface for Waveform reco algorithms
class IWaveformReconstructionTool : virtual public IAlgTool
@@ -31,6 +32,11 @@ class IWaveformReconstructionTool : virtual public IAlgTool
virtual ~IWaveformReconstructionTool() = default;
+ // Reconstruct all waveforms
+ virtual StatusCode reconstructAll(const RawWaveformContainer& waveContainer,
+ const xAOD::WaveformClock* clock,
+ xAOD::WaveformHitContainer* container) const = 0;
+
// Reconstruct all peaks in a raw waveform
virtual StatusCode reconstruct(const RawWaveform& wave,
const xAOD::WaveformClock* clock,
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
index 5de393166579eef162e308acd31b0b12b6a12039..69a1d2be1f1b4a9dae3f230bf9c9a4915d3b18fe 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
@@ -10,8 +10,6 @@
#include "WaveformReconstructionTool.h"
-#include "xAODFaserWaveform/WaveformHit.h"
-
#include "TH1F.h"
#include "TF1.h"
#include "TFitResult.h"
@@ -42,6 +40,215 @@ WaveformReconstructionTool::initialize() {
}
// Reconstruction step
+StatusCode
+WaveformReconstructionTool::reconstructAll(
+ const RawWaveformContainer& waveContainer,
+ const xAOD::WaveformClock* clock,
+ xAOD::WaveformHitContainer* hitContainer) const {
+
+ ATH_MSG_DEBUG(" reconstructAll called ");
+
+ // Reconstruct each waveform
+ for( const auto& wave : waveContainer) {
+
+ ATH_MSG_DEBUG("Reconstruct waveform for channel " << wave->channel());
+
+ // Reconstruct the hits, may be more than one, so pass container
+ CHECK( this->reconstruct(*wave, clock, hitContainer) );
+ }
+
+ if (m_ensureChannelHits) {
+ ATH_MSG_DEBUG("Ensure all channels have hits at peak time");
+ ensureHits(waveContainer, clock, hitContainer);
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+//
+// Make sure we have a hit for each channel at the time when
+// there is a significant pulse found in the detector
+//
+void
+WaveformReconstructionTool::ensureHits(
+ const RawWaveformContainer& waveContainer,
+ const xAOD::WaveformClock* clock,
+ xAOD::WaveformHitContainer* hitContainer) const {
+
+ ATH_MSG_DEBUG(" ensureHits called ");
+
+ // Find peak time (most significant hit)
+ xAOD::WaveformHit* peakHit = NULL;
+
+ for( const auto& hit : *hitContainer) {
+
+ if (peakHit == NULL) {
+ peakHit = hit;
+ } else {
+ if ( hit->peak() > peakHit->peak() ) peakHit = hit;
+ }
+
+ }
+
+ // Didn't find anything?
+ if (peakHit == NULL) return;
+ if (peakHit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) return;
+
+ ATH_MSG_DEBUG("Found peak hit in channel " << peakHit->channel() << " at time " << peakHit->localtime());
+
+ // Now go through all of the channels and check if there is a hit
+ // close in time to the peakHit
+ for( const auto& wave : waveContainer) {
+
+ // Don't worry about the peak channel, we know this has a hit...
+ if (wave->channel() == peakHit->channel()) continue;
+
+ ATH_MSG_DEBUG("Checking for hit in channel " << wave->channel());
+
+ bool found = false;
+ // Look for a baseline-only hit that we can update
+ xAOD::WaveformHit* baselineHit = NULL;
+
+ // There aren't so many hits, just loop over container
+ for( const auto& hit : *hitContainer) {
+ if (hit->channel() != wave->channel()) continue;
+
+ // Is this above threshold?
+ if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) {
+ baselineHit = hit;
+ continue;
+ }
+
+ // OK, this is the right channel, check the time
+ float dtime = abs(hit->localtime() - peakHit->localtime());
+ if (dtime > m_hitTimeDifference) continue;
+
+ // We have found a hit in the right channel at the right time
+ found = true;
+ ATH_MSG_DEBUG("Found hit in channel " << hit->channel()
+ << " at time " << hit->localtime());
+ break;
+ }
+
+ // Is there a hit? If so, go to next waveform/channel
+ if (found) continue;
+
+ ATH_MSG_DEBUG("No hit found for channel " << wave->channel()
+ << " at time " << peakHit->localtime());
+
+ // Do we have a baseline-only hit we can use?
+ xAOD::WaveformHit* newhit = NULL;
+ if (baselineHit == NULL) {
+ // No, make a new hit here
+ newhit = new xAOD::WaveformHit();
+ hitContainer->push_back(newhit);
+
+ // Mark this as a secondary hit
+ newhit->set_status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED);
+ newhit->set_status_bit(xAOD::WaveformStatus::SECONDARY);
+
+ // Set digitizer channel and identifier
+ newhit->set_channel(wave->channel());
+ newhit->set_id(wave->identify32().get_compact());
+
+ // Make sure we have ADC counts
+ if (wave->adc_counts().size() == 0) {
+ ATH_MSG_WARNING( "Found waveform for channel " << wave->channel()
+ << " with size " << wave->adc_counts().size() << "!");
+
+ newhit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_MISSING);
+ continue;
+ }
+
+ if (wave->adc_counts().size() != wave->n_samples()) {
+ ATH_MSG_WARNING( "Found waveform for channel " << wave->channel()
+ << " with size " << wave->adc_counts().size()
+ << " not equal to number of samples " << wave->n_samples());
+
+ newhit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_INVALID);
+ continue;
+ }
+
+ findBaseline(*wave, newhit);
+
+ } else {
+ // Use the existing baseline hit
+ newhit = baselineHit;
+ }
+
+ // Check for problems
+ if (newhit->status_bit(xAOD::WaveformStatus::BASELINE_FAILED)) continue;
+
+ // Set range for windowed data
+ unsigned int lo_edge = peakHit->time_vector().front()/2.;
+ unsigned int hi_edge = peakHit->time_vector().back()/2.;
+
+ ATH_MSG_DEBUG("Windowing waveform from " << lo_edge << " to " << hi_edge);
+ std::vector<float> wtime(hi_edge-lo_edge+1);
+ std::vector<float> wwave(hi_edge-lo_edge+1);
+ for (unsigned int i=lo_edge; i<=hi_edge; i++) {
+ unsigned int j = i-lo_edge;
+ wtime[j] = 2.*i;
+ wwave[j] = newhit->baseline_mean() - wave->adc_counts()[i];
+ //ATH_MSG_DEBUG(" Time: " << wtime[j] << " Wave: " << wwave[j]);
+ }
+
+ newhit->set_time_vector(wtime);
+ newhit->set_wave_vector(wwave);
+
+ //
+ // Find some raw values
+ WaveformFitResult raw = findRawHitValues(wtime, wwave);
+ newhit->set_peak(raw.peak);
+ newhit->set_mean(raw.mean);
+ newhit->set_width(raw.sigma);
+ newhit->set_integral(raw.integral);
+ newhit->set_localtime(raw.mean);
+ newhit->set_raw_peak(raw.peak);
+ newhit->set_raw_integral(raw.integral);
+
+ //
+ // Find time from clock
+ if (!clock || (clock->frequency() <= 0.)) {
+ newhit->set_status_bit(xAOD::WaveformStatus::CLOCK_INVALID);
+ newhit->set_bcid_time(-1.);
+ } else {
+ newhit->set_bcid_time(clock->time_from_clock(newhit->localtime()));
+ }
+
+ } // End of loop over waveContainer
+}
+
+// Find the baseline
+WaveformBaselineData&
+WaveformReconstructionTool::findBaseline(const RawWaveform& raw_wave,
+ xAOD::WaveformHit* hit) const {
+
+ //
+ // Find baseline
+ static WaveformBaselineData baseline;
+
+ if (m_useSimpleBaseline.value())
+ baseline = findSimpleBaseline(raw_wave);
+ else
+ baseline = findAdvancedBaseline(raw_wave);
+
+ if (!(baseline.valid)) {
+ ATH_MSG_WARNING("Failed to reconstruct baseline!");
+
+ hit->set_baseline_mean(0.);
+ hit->set_baseline_rms(-1.);
+ hit->set_status_bit(xAOD::WaveformStatus::BASELINE_FAILED);
+
+ } else {
+ // Save baseline to hit collection object
+ hit->set_baseline_mean(baseline.mean);
+ hit->set_baseline_rms(baseline.rms);
+ }
+
+ return baseline;
+}
+
StatusCode
WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave,
const xAOD::WaveformClock* clock,
@@ -64,16 +271,15 @@ WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave,
hit->set_channel(raw_wave.channel());
hit->set_id(raw_wave.identify32().get_compact());
- //
- // Make some sanity checks on the waveform data
+ // Make sure we have ADC counts
if (raw_wave.adc_counts().size() == 0) {
ATH_MSG_WARNING( "Found waveform for channel " << raw_wave.channel()
<< " with size " << raw_wave.adc_counts().size() << "!");
hit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_MISSING);
return StatusCode::SUCCESS;
- }
-
+ }
+
if (raw_wave.adc_counts().size() != raw_wave.n_samples()) {
ATH_MSG_WARNING( "Found waveform for channel " << raw_wave.channel()
<< " with size " << raw_wave.adc_counts().size()
@@ -83,28 +289,12 @@ WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave,
return StatusCode::SUCCESS;
}
- //
- // Find baseline
- WaveformBaselineData baseline;
-
- if (m_useSimpleBaseline.value())
- baseline = findSimpleBaseline(raw_wave);
- else
- baseline = findAdvancedBaseline(raw_wave);
-
- if (!(baseline.valid)) {
- ATH_MSG_WARNING("Failed to reconstruct baseline!");
+ // Find the baseline
+ WaveformBaselineData baseline = findBaseline(raw_wave, hit);
- hit->set_baseline_mean(0.);
- hit->set_baseline_rms(-1.);
- hit->set_status_bit(xAOD::WaveformStatus::BASELINE_FAILED);
-
- return StatusCode::SUCCESS;
- }
-
- // Save baseline to hit collection object
- hit->set_baseline_mean(baseline.mean);
- hit->set_baseline_rms(baseline.rms);
+ // Check that we have data to work with
+ // If any status bits are set, this is bad
+ if (hit->status()) return StatusCode::SUCCESS;
//
// Create baseline-subtracted data array for both time and signal
@@ -349,10 +539,19 @@ WaveformReconstructionTool::findAdvancedBaseline(const RawWaveform& raw_wave) co
TH1F h1("", "", nbins, xlo, xhi);
// Fill this histogram with the waveform
+ // Only use a subset, based on where we think there won't be signal
+ const std::vector<unsigned int>& counts = raw_wave.adc_counts();
+ for (int i=m_baselineSampleLo; i<=m_baselineSampleHi; i++) {
+ //ATH_MSG_INFO( "Entry " << i << " Value " << counts[i] );
+ h1.Fill(counts[i]);
+ }
+
+ /*
for (auto value : raw_wave.adc_counts()) {
//ATH_MSG_DEBUG( "Found value " << value );
h1.Fill(value);
}
+ */
// Find max bin
int maxbin = h1.GetMaximumBin();
@@ -366,8 +565,8 @@ WaveformReconstructionTool::findAdvancedBaseline(const RawWaveform& raw_wave) co
ATH_MSG_DEBUG( "Filling 2nd histogram from " << xlo << " to " << xhi);
TH1F h2("", "", nbins, xlo, xhi);
- for (auto value : raw_wave.adc_counts()) {
- h2.Fill(value);
+ for (int i=m_baselineSampleLo; i<=m_baselineSampleHi; i++) {
+ h2.Fill(counts[i]);
}
// Start with full histogram range
@@ -459,7 +658,7 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons
rfit.integral = 2*tot; // Factor of 2 because at 500 MHz, dt = 2 ns
rfit.time = rfit.mean;
- ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak);
+ ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak << " Integral: " << rfit.integral);
// Fix bad values
if (isnan(rfit.sigma)) {
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.h b/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
index 1bb0f20d4d2901468fa9982fd7bbed26c9c2d73b..e7e61205f62af237aacaa2d7b8f803638856d7b0 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
@@ -14,6 +14,8 @@
#include "WaveRecTools/IWaveformReconstructionTool.h"
#include "WaveRawEvent/RawWaveform.h"
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "xAODFaserWaveform/WaveformHit.h"
#include "WaveformBaselineData.h"
#include "WaveformFitResult.h"
@@ -35,25 +37,40 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
/// Retrieve the necessary services in initialize
StatusCode initialize();
+ /// Reconstruct all hits from waveform container
+ virtual StatusCode reconstructAll(const RawWaveformContainer& waveContainer,
+ const xAOD::WaveformClock* clock,
+ xAOD::WaveformHitContainer* hitContainer) const;
+
/// Reconstruct hits from waveform
virtual StatusCode reconstruct(const RawWaveform& wave,
const xAOD::WaveformClock* clock,
- xAOD::WaveformHitContainer* container) const;
+ xAOD::WaveformHitContainer* hitContainer) const;
private:
//
// Baseline Estimation Parameters
BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", false};
+
+ // Minimum number of samples needed to calculate simple baseline
+ // Just average these first n values
IntegerProperty m_samplesForBaselineAverage{this, "SamplesForBaselineAverage", 40};
+
//
- // Parameters of initial histogram to find baseline max
+ // Parameters for advanced baseline determination
+ // Parameters of initial histogram to find most likely ADC reading
// Range and bins to use, ratio should be an integer (bin width)
IntegerProperty m_baselineRange{this, "BaselineRange", 16000};
IntegerProperty m_baselineRangeBins{this, "BaselineRangeBins", 320};
+ // Range of samples to use to find baseline (to avoid signal contamination)
+ // In samples (2 ns intervals) so 350 -> 700 ns in digitizer window
+ IntegerProperty m_baselineSampleLo{this, "BaselineSampleLo", 0};
+ IntegerProperty m_baselineSampleHi{this, "BaselineSampleHi", 350};
+
//
// Parameters for the Gaussian fit to the baseline peak (in counts)
// Range is total range to use to find truncated mean and width
@@ -82,9 +99,19 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
// Fraction of peak to set local hit time
FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.45};
+ //
+ // Ensure each channel has a waveform hit at time of most significant
+ // hit in the event
+ BooleanProperty m_ensureChannelHits{this, "EnsureChannelHits", true};
+ // Max Time difference in ns to say a hit exists in a different channel
+ FloatProperty m_hitTimeDifference{this, "HitTimeDifference", 10.};
+
// Baseline algorithms
WaveformBaselineData& findSimpleBaseline(const RawWaveform& wave) const;
WaveformBaselineData& findAdvancedBaseline(const RawWaveform& wave) const;
+ WaveformBaselineData& findBaseline(const RawWaveform& wave,
+ xAOD::WaveformHit* hit) const;
+
// Find peak in wave, return windowed region in windowed_time and windowed_wave
// Windowed region is removed from original vectors
@@ -112,6 +139,11 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
const std::vector<float> wave) const;
+ /// Create hit in all channels at time of peak signal
+ void ensureHits(const RawWaveformContainer& waveContainer,
+ const xAOD::WaveformClock* clock,
+ xAOD::WaveformHitContainer* hitContainer) const;
+
};
#endif // WAVERECTOOLS_WAVEFORMRECONSTRUCTIONTOOL_H
diff --git a/faser-common b/faser-common
index 89ce6a07128eb2ebc367b6b68f29c9c88220e3e6..e7a3a145593c7bbda73b64a7106ffe3e43f8ff94 160000
--- a/faser-common
+++ b/faser-common
@@ -1 +1 @@
-Subproject commit 89ce6a07128eb2ebc367b6b68f29c9c88220e3e6
+Subproject commit e7a3a145593c7bbda73b64a7106ffe3e43f8ff94