diff --git a/Database/ConnectionManagement/FaserAuthentication/data/dblookup.xml b/Database/ConnectionManagement/FaserAuthentication/data/dblookup.xml
index 52542a2837125d540824617b1085a5d710c3507f..9cebc2e82b233d6191d7671cf6d9d6997cad23f9 100644
--- a/Database/ConnectionManagement/FaserAuthentication/data/dblookup.xml
+++ b/Database/ConnectionManagement/FaserAuthentication/data/dblookup.xml
@@ -26,4 +26,9 @@
    <service name="sqlite_file:///cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/sqlite200/CABP200.db" accessMode="read" />
 </logicalservice>
 
+<logicalservice name="COOLOFL_TRIGGER">
+   <service name="sqlite_file:data/sqlite200/waveform_reco.db" accessMode="read" />
+   <service name="sqlite_file:///cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/sqlite200/waveform_reco.db" accessMode="read" />
+</logicalservice>
+
 </servicelist>
diff --git a/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h b/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h
index 9581de5530a2c65c8ddddeb870c971a14b988b34..a1f42f6e86db846ea544daedb482ba8aba8fbdf5 100644
--- a/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h
+++ b/Waveform/WaveRawEvent/WaveRawEvent/RawWaveform.h
@@ -62,6 +62,7 @@ public:
   // Waveform data
   unsigned int channel() const;
   const std::vector<unsigned int>& adc_counts() const;
+  size_t size() const {return m_adc_counts.size();}
 
   // Return channel identifier
   Identifier identify() const;
diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
index fb1cc4d4261adc110578d959a9fa3c24801c8b77..ad44210937a7f4aad1859f3c7a8100d55d4070c3 100644
--- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
+++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py
@@ -6,6 +6,7 @@ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
 from AthenaConfiguration.ComponentFactory import CompFactory
 
 from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+from WaveformConditionsTools.WaveformTimingConfig import WaveformTimingCfg
 
 WaveformReconstructionTool = CompFactory.WaveformReconstructionTool
 ClockReconstructionTool = CompFactory.ClockReconstructionTool
@@ -32,6 +33,8 @@ def WaveformReconstructionCfg(flags, naive = False):
     acc.merge(WaveformHitRecCfg(flags, "PreshowerWaveformRecAlg", "Preshower"))
     acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo"))
 
+    acc.merge(WaveformTimingCfg(flags))
+
     return acc
 
 # Return configured WaveformClock reconstruction algorithm
diff --git a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
index 3248939a0879178789e7c770fb584fc0b85cce07..243b88189350fb12bf5aab958a4d1ba7956a6a88 100644
--- a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
+++ b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
@@ -31,9 +31,10 @@ RawWaveformRecAlg::finalize() {
 
   ATH_MSG_INFO( m_numberOfEvents << " events processed" );
   if ( m_numberOfEvents > 0) { 
-    ATH_MSG_INFO( m_numberOfWaveforms << " waveforms found" );
-    ATH_MSG_INFO( m_numberOfOverflows << " overflows" );
-    ATH_MSG_INFO( m_numberOfFitErrors << " fit errors" );
+    ATH_MSG_INFO( m_numberOfWaveforms   << " waveforms found over threshold" );
+    ATH_MSG_INFO( m_numberOfSecondaries << " secondary waveforms found" );
+    ATH_MSG_INFO( m_numberOfOverflows   << " overflows" );
+    ATH_MSG_INFO( m_numberOfFitErrors   << " fit errors" );
   }
 
   return StatusCode::SUCCESS;
@@ -64,30 +65,45 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
     return StatusCode::SUCCESS;
   }
 
+  // First reconstruct the primary hit (based on trigger time)
+  for( const auto& wave : *waveformHandle) {
+    ATH_MSG_DEBUG("Reconstruct primary waveform for channel " << wave->channel());
+    CHECK( m_recoTool->reconstructPrimary(*wave, hitContainerHandle.ptr()) );
+  }
+
+  // Second, reconstruct any additional out of time hits
+  if (m_findMultipleHits) {
+    for( const auto& wave : *waveformHandle) {
+      ATH_MSG_DEBUG("Reconstruct secondary waveform for channel " << wave->channel());
+      CHECK( m_recoTool->reconstructSecondary(*wave, hitContainerHandle.ptr()) );
+    }
+  }
+
   // Also find the clock information
   SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_clockKey, ctx);
   const xAOD::WaveformClock* clockptr = NULL;
 
+  // Fix timing for all hits
   // Can survive without this, but make a note
   if ( clockHandle.isValid() ) {
     ATH_MSG_DEBUG("Found ReadHandle for WaveformClock");
     clockptr = clockHandle.ptr();
+    CHECK( m_recoTool->setLocalTime(clockptr, hitContainerHandle.ptr()) );
   } else {
     ATH_MSG_WARNING("Didn't find ReadHandle for WaveformClock!");
   }
 
-
-  // Reconstruct all waveforms
-  CHECK( m_recoTool->reconstructAll(*waveformHandle, clockptr, hitContainerHandle.ptr()) );
-
   ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' filled with "<< hitContainerHandle->size() <<" items");
 
   // Keep track of some statistics
   m_numberOfEvents++;
   for (const auto& hit : *(hitContainerHandle.ptr())) {
     if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) continue;
+
     m_numberOfWaveforms++;
     if (hit->status_bit(xAOD::WaveformStatus::WAVE_OVERFLOW)) m_numberOfOverflows++;
+    if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) m_numberOfSecondaries++;
+    
     if (hit->status_bit(xAOD::WaveformStatus::GFIT_FAILED)) {
       m_numberOfFitErrors++;
     } else if (hit->status_bit(xAOD::WaveformStatus::CBFIT_FAILED)) {
diff --git a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.h b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.h
index e57501a730dd5f67c38435cb9594b0a7def2ec75..4de415d21d162a3aa2b64fe7f1ef1bb14e0237d9 100644
--- a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.h
+++ b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.h
@@ -42,6 +42,10 @@ class RawWaveformRecAlg : public AthReentrantAlgorithm {
   virtual StatusCode finalize() override;
   //@}
 
+  //
+  // Look for more than one hit in each channel
+  BooleanProperty m_findMultipleHits{this, "FindMultipleHits", true};
+
  private:
 
   /** @name Disallow default instantiation, copy, assignment */
@@ -90,6 +94,7 @@ class RawWaveformRecAlg : public AthReentrantAlgorithm {
   //@{
   mutable std::atomic<int> m_numberOfEvents{0};
   mutable std::atomic<int> m_numberOfWaveforms{0};
+  mutable std::atomic<int> m_numberOfSecondaries{0};
   mutable std::atomic<int> m_numberOfOverflows{0};
   mutable std::atomic<int> m_numberOfFitErrors{0};
   //@}
diff --git a/Waveform/WaveRecTools/CMakeLists.txt b/Waveform/WaveRecTools/CMakeLists.txt
index d8f3e6f053232477a1b1501a26e3a88efe9272c1..f7f9672886ef12ab2fd9441447bf0f0a0f583418 100644
--- a/Waveform/WaveRecTools/CMakeLists.txt
+++ b/Waveform/WaveRecTools/CMakeLists.txt
@@ -13,13 +13,16 @@ atlas_add_library( WaveRecToolsLib
                    WaveRecTools/*.h src/*.cxx src/*.h
                    PUBLIC_HEADERS WaveRecTools
                    PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
-                   LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent xAODFaserWaveform
+                   LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives 
+		   WaveformConditionsToolsLib WaveRawEvent xAODFaserWaveform
                    PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
 		   )
 
 atlas_add_component( WaveRecTools
 		     src/components/*.cxx 
 		     INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
-                     LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveRecToolsLib )
+                     LINK_LIBRARIES ${ROOT_LIBRARIES} 
+		     WaveformConditionsToolsLib AthenaBaseComps GaudiKernel 
+		     WaveRecToolsLib)
 
 
diff --git a/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h b/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
index cc10197b262f0695327f7304a84998d198300678..c72c2502c9391861a3aec0696f47d568a7d7aa2b 100644
--- a/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
+++ b/Waveform/WaveRecTools/WaveRecTools/IWaveformReconstructionTool.h
@@ -32,15 +32,17 @@ 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, 
-				 xAOD::WaveformHitContainer* container) const = 0;
+  // Reconstruct hits in trigger window
+  virtual StatusCode reconstructPrimary(const RawWaveform& wave, 
+					xAOD::WaveformHitContainer* container) const = 0;
+
+  // Reconstruct secondary hits anywhere in the waveform
+  virtual StatusCode reconstructSecondary(const RawWaveform& wave, 
+					  xAOD::WaveformHitContainer* container) const = 0;
+
+  // Set local hit times from LHC clock
+  virtual StatusCode setLocalTime(const xAOD::WaveformClock* clock,
+				  xAOD::WaveformHitContainer* container) const = 0;
 
 };
 
diff --git a/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx b/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
index c45ea51007362eecd0ab7827f50123c05ddebab7..825cc755862200c0dd008eaf3c273cb23b37f83a 100644
--- a/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/ClockReconstructionTool.cxx
@@ -102,10 +102,10 @@ ClockReconstructionTool::reconstruct(const RawWaveform& raw_wave,
       ATH_MSG_DEBUG("Index: " << i << " Freq: " << i*freqmult << " Mag: " << magnitude[i]);
   }
 
-  // Store results
-  clockdata->set_dc_offset(magnitude[0]);
+  // Store results (amplitides in mV)
+  clockdata->set_dc_offset(raw_wave.mv_per_bit()*magnitude[0]);
+  clockdata->set_amplitude(raw_wave.mv_per_bit()*magnitude[imax]);
   clockdata->set_frequency(imax * freqmult);
-  clockdata->set_amplitude(magnitude[imax]);
   clockdata->set_phase(atan2(im_full[imax], re_full[imax])); // Not a bug, atan2(y,x)!
 
   ATH_MSG_DEBUG("Before correcting for finite resolution:");
@@ -133,7 +133,7 @@ ClockReconstructionTool::reconstruct(const RawWaveform& raw_wave,
 
   clockdata->set_frequency( (imax+dm) * freqmult );
   clockdata->set_phase (phase);
-  clockdata->set_amplitude( 2*M_PI*dm*magnitude[imax] / sin(M_PI * dm) );
+  clockdata->set_amplitude( raw_wave.mv_per_bit() * 2*M_PI*dm*magnitude[imax] / sin(M_PI * dm) );
 
   ATH_MSG_DEBUG("After correcting for finite resolution:");
   ATH_MSG_DEBUG(*clockdata);
diff --git a/Waveform/WaveRecTools/src/ClockReconstructionTool.h b/Waveform/WaveRecTools/src/ClockReconstructionTool.h
index ea7ec2ec7e2a28e49d587f9a7f072678ee4bdbc5..78e6ea1770c069314ee15a068458f88623b63bda 100644
--- a/Waveform/WaveRecTools/src/ClockReconstructionTool.h
+++ b/Waveform/WaveRecTools/src/ClockReconstructionTool.h
@@ -51,8 +51,8 @@ class ClockReconstructionTool: public extends<AthAlgTool, IClockReconstructionTo
   void checkResult(const RawWaveform& raw_wave,
 		   xAOD::WaveformClock* clockdata) const;
 
-  // Limits to print warnings
-  FloatProperty m_amplitude_min{this, "AmplitudeMin", 1000.};
+  // Limits to print warnings (amplitude in mV)
+  FloatProperty m_amplitude_min{this, "AmplitudeMin", 500.};
   FloatProperty m_frequency_min{this, "FrequencyMin", 40.0};
   FloatProperty m_frequency_max{this, "FrequencyMax", 40.1};
 
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
index 04f857e038fa4c84ca2cefd88b58a0f153758752..296b36d2090c569194ddcecde8a3ba5cdbce2ec0 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.cxx
@@ -39,188 +39,277 @@ WaveformReconstructionTool::initialize() {
   } else {
     ATH_MSG_INFO("Will use fit to determine baseline");
   }
+
+  ATH_CHECK( m_timingTool.retrieve() );
+
   return StatusCode::SUCCESS;
 }
 
-// Reconstruction step
+//
+// Form primary hits using trigger time
+//
 StatusCode
-WaveformReconstructionTool::reconstructAll(
-	 const RawWaveformContainer& waveContainer, 
-	 const xAOD::WaveformClock* clock, 
+WaveformReconstructionTool::reconstructPrimary(
+	 const RawWaveform& wave,
 	 xAOD::WaveformHitContainer* hitContainer) const {
 
-  ATH_MSG_DEBUG(" reconstructAll called ");
+  ATH_MSG_DEBUG(" reconstructPrimary called");
 
-  // Reconstruct each waveform
-  for( const auto& wave : waveContainer) {
+  // Get the nominal trigger time (in ns) from config
+  float trigger_time = m_timingTool->nominalTriggerTime();
 
-    ATH_MSG_DEBUG("Reconstruct waveform for channel " << wave->channel());
+  xAOD::WaveformHit* newhit = new xAOD::WaveformHit();
+  hitContainer->push_back(newhit);
 
-    // Reconstruct the hits, may be more than one, so pass container
-    CHECK( this->reconstruct(*wave, clock, hitContainer) );
+  // 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);
+    return StatusCode::SUCCESS;
+  } 
+
+  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);
+    return StatusCode::SUCCESS;
   }
 
-  if (m_ensureChannelHits) {
-    ATH_MSG_DEBUG("Ensure all channels have hits at peak time");
-    ensureHits(waveContainer, clock, hitContainer);
+  // Find the baseline for this waveform
+  findBaseline(wave, newhit);
+
+  // Check for problems
+  if (newhit->status_bit(xAOD::WaveformStatus::BASELINE_FAILED)) 
+    return StatusCode::SUCCESS;
+
+  // Set range for windowed data in digitizer samples
+  float offset = m_timingTool->triggerTimeOffset(wave.channel());
+
+  int lo_edge = int((trigger_time+offset)/2.) + m_windowStart;
+  int hi_edge = int((trigger_time+offset)/2.) + m_windowStart + m_windowWidth;
+
+  // Fill raw hit values
+  fillRawHitValues(wave, lo_edge, hi_edge, newhit);
+
+  // Check if this is over threshold
+  if (newhit->peak() < newhit->baseline_rms() * m_primaryPeakThreshold) {
+    ATH_MSG_DEBUG("Primary hit failed threshold");
+    newhit->set_status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED);
+  } else {
+    // Reconstruct hit in this range
+    reconstructHit(newhit);
   }
 
   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
+// Form primary hits using trigger time
 //
-void
-WaveformReconstructionTool::ensureHits(
-	 const RawWaveformContainer& waveContainer, 
-	 const xAOD::WaveformClock* clock, 
+StatusCode
+WaveformReconstructionTool::reconstructSecondary(
+	 const RawWaveform& wave,
 	 xAOD::WaveformHitContainer* hitContainer) const {
 
-  ATH_MSG_DEBUG(" ensureHits called ");
+  ATH_MSG_DEBUG(" reconstructSecondary called");
 
-  // Find peak time (most significant hit)
-  xAOD::WaveformHit* peakHit = NULL;
+  // Find existing hit for this channel to get baseline
+  xAOD::WaveformHit* primaryHit = NULL;
 
   for( const auto& hit : *hitContainer) {
 
-    if (peakHit == NULL) {
-      peakHit = hit;
-    } else {
-      if ( hit->peak() > peakHit->peak() ) peakHit = hit;
+    if (hit->channel() == wave.channel()) {
+      primaryHit = hit;
+      break;
     }
+  }
 
+  // Did we find the primary hit for this channel?
+  if (!primaryHit) {
+    ATH_MSG_ERROR("found no primary hit for channel " << wave.channel() << "!");
+    return StatusCode::FAILURE;
   }
 
-  // Didn't find anything?
-  if (peakHit == NULL) return;
-  if (peakHit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) return;
+  WaveformBaselineData baseline;
 
-  ATH_MSG_DEBUG("Found peak hit in channel " << peakHit->channel() << " at time " << peakHit->localtime());
+  baseline.mean = primaryHit->baseline_mean();
+  baseline.rms = primaryHit->baseline_rms();
+  
+  // Find the secondary peak position
+  int ipeak;
 
-  // 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) {
+  // Is there already a peak in the primary?
+  if (primaryHit->threshold()) {
 
-    // Don't worry about the peak channel, we know this has a hit...
-    if (wave->channel() == peakHit->channel()) continue;
+    ATH_MSG_DEBUG("Looking for secondary hit with primary hit above threshold");
 
-    ATH_MSG_DEBUG("Checking for hit in channel " << wave->channel());
+    // Look before and after window
+    int lo_edge = int(primaryHit->time_vector().front()/2.);
+    int hi_edge = int(primaryHit->time_vector().back()/2.);
 
-    bool found = false;
-    // Look for a baseline-only hit that we can update
-    xAOD::WaveformHit* baselineHit = NULL;
+    std::vector<float> wwave_lo(lo_edge);
+    std::vector<float> wwave_hi(wave.adc_counts().size() - hi_edge - 1);
 
-    // There aren't so many hits, just loop over container
-    for( const auto& hit : *hitContainer) {
-      if (hit->channel() != wave->channel()) continue;
+    int ipeak_lo = -1.;
+    int ipeak_hi = -1.;
 
-      // Is this above threshold?
-      if (hit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) {
-	baselineHit = hit;
-	continue;
+    // Look before
+    if (m_findSecondaryBefore) {
+      for (int i=0; i<lo_edge; i++) {
+	wwave_lo[i] = baseline.mean - wave.mv_per_bit() * wave.adc_counts()[i];
       }
 
-      // OK, this is the right channel, check the time
-      float dtime = abs(hit->localtime() - peakHit->localtime());
-      if (dtime > m_hitTimeDifference) continue;
+      ipeak_lo = findPeak(baseline, m_secondaryPeakThreshold, wwave_lo);
 
-      // 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;
+      if (ipeak_lo < 0) {
+	ATH_MSG_DEBUG("No hit found before " << lo_edge);
+      } else {
+	ATH_MSG_DEBUG("Hit found at " << ipeak_lo << " before " << lo_edge);
+      }
     }
 
-    // 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);
+    // Look after
+    if (m_findSecondaryAfter) {
+      for (unsigned int i=(hi_edge+1); i<wave.adc_counts().size(); i++) {
+	wwave_hi[(i-(hi_edge+1))] = baseline.mean - wave.mv_per_bit() * wave.adc_counts()[i];
+      }
 
-      // Mark this as a secondary hit
-      newhit->set_status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED);
-      newhit->set_status_bit(xAOD::WaveformStatus::SECONDARY);
+      ipeak_hi = findPeak(baseline, m_secondaryPeakThreshold, wwave_hi);
 
-      // Set digitizer channel and identifier
-      newhit->set_channel(wave->channel());
-      newhit->set_id(wave->identify32().get_compact());
+      // Is this too close to the primary hit?
+      if (ipeak_hi < 5) {
+	ATH_MSG_DEBUG("Found hit after at " << (ipeak_hi + hi_edge + 1)<< " but too close to edge");
+	ipeak_hi = -1;
+      }
 
-      // 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;
+      if (ipeak_hi < 0) {
+	ATH_MSG_DEBUG("No hit found after " << hi_edge);
+      } else {
+	ATH_MSG_DEBUG("Hit found at " << ipeak_hi << " after " << hi_edge);
       }
+    }
 
-      findBaseline(*wave, newhit);
+    // Nothing found
+    if (ipeak_lo < 0 && ipeak_hi < 0)
+      return StatusCode::SUCCESS;
+
+    // Both?
+    if (ipeak_lo >= 0 && ipeak_hi >= 0) {
+
+      // Pick the largest signal
+      if (wwave_lo[ipeak_lo] >= wwave_hi[ipeak_hi]) {
+	ipeak = ipeak_lo;
+	ATH_MSG_DEBUG("Picked before as " << wwave_lo[ipeak_lo]
+		      << " > " << wwave_hi[ipeak_hi]);
+      } else {
+	ipeak = ipeak_hi + hi_edge + 1;
+	ATH_MSG_DEBUG("Picked after as " << wwave_lo[ipeak_lo]
+		      << " < " << wwave_hi[ipeak_hi]);
+      }
 
+    
+    } else if (ipeak_lo > 0) {
+      ipeak = ipeak_lo;
+      ATH_MSG_DEBUG("Peak before with " << wwave_lo[ipeak_lo]);
     } else {
-      // Use the existing baseline hit
-      newhit = baselineHit;
+      ATH_MSG_DEBUG("Peak after with " << wwave_hi[ipeak_hi]);
+      ipeak = ipeak_hi+hi_edge+1;
     }
 
-    // Check for problems
-    if (newhit->status_bit(xAOD::WaveformStatus::BASELINE_FAILED)) continue;
+  } else {
 
-    // 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("Looking for secondary hit without primary hit above threshold");
+    std::vector<float> wwave(wave.adc_counts().size());
+    for (unsigned int i=0; i<wave.adc_counts().size(); i++) {
+      wwave[i] = baseline.mean - wave.mv_per_bit() * wave.adc_counts()[i];
+    }
 
-    ATH_MSG_DEBUG("Windowing waveform from " << lo_edge << " to " << hi_edge);
+    ipeak = findPeak(baseline, m_secondaryPeakThreshold, wwave);
 
-    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; // 2 ns per sample at 500 MHz sampling
-      wwave[j] = newhit->baseline_mean() - wave->mv_per_bit()*wave->adc_counts()[i];
-      //ATH_MSG_DEBUG(" Time: " << wtime[j] << " Wave: " << wwave[j]);
-    }
+    // Nothing found
+    if (ipeak < 0) 
+      return StatusCode::SUCCESS;
+
+    ATH_MSG_DEBUG("Found secondary peak with no primary " << wwave[ipeak]);
+  }
 
-    newhit->set_time_vector(wtime);
-    newhit->set_wave_vector(wwave);
+  // We seem to have a secondary hit
+  xAOD::WaveformHit* newhit = new xAOD::WaveformHit();
+  hitContainer->push_back(newhit);
 
-    //
-    // 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);
+  // Fill values
+  newhit->set_channel(wave.channel());
+  newhit->set_id(wave.identify32().get_compact());
+  newhit->set_status_bit(xAOD::WaveformStatus::SECONDARY);
+  newhit->set_baseline_mean(baseline.mean);
+  newhit->set_baseline_rms(baseline.rms);
+
+  // Set range for windowed data in digitizer samples
+  int lo_edge = ipeak + m_windowStart;
+  int hi_edge = ipeak + m_windowStart + m_windowWidth;
+
+  // Fill raw hit values
+  fillRawHitValues(wave, lo_edge, hi_edge, newhit);
+
+  // Must be over threshold, so reconstruct here
+  reconstructHit(newhit);
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode
+WaveformReconstructionTool::setLocalTime(const xAOD::WaveformClock* clock,
+					 xAOD::WaveformHitContainer* container) const {
+
+  ATH_MSG_DEBUG(" setLocalTime called ");
+
+  // Check the container
+  if (!container) {
+    ATH_MSG_ERROR("WaveformHitCollection passed to setLocalTime() is null!");
+    return StatusCode::FAILURE;
+  }
+
+  bool clock_valid;
+
+  //
+  // Find time from clock
+  if (!clock || (clock->frequency() <= 0.)) {
+    clock_valid = false;
+  } else {
+    clock_valid = true;
+  }
+
+  float trigger_time = m_timingTool->nominalTriggerTime();
+  float offset;
+
+  // Should actually find the time of the trigger here 
+  // and set bcid time offset from that
+  // Loop through hits and set local time
+  for( const auto& hit : *container) {
 
     //
     // Find time from clock
-    if (!clock || (clock->frequency() <= 0.)) {
-      newhit->set_status_bit(xAOD::WaveformStatus::CLOCK_INVALID);
-      newhit->set_bcid_time(-1.);
+    if (clock_valid) {
+      hit->set_bcid_time(clock->time_from_clock(hit->localtime()));
     } else {
-      newhit->set_bcid_time(clock->time_from_clock(newhit->localtime()));
+      hit->set_status_bit(xAOD::WaveformStatus::CLOCK_INVALID);
+      hit->set_bcid_time(-1.);
     }
 
-  } // End of loop over waveContainer
+    // Also set time with respect to nominal trigger
+    offset = m_timingTool->triggerTimeOffset(hit->channel());
+    hit->set_trigger_time(hit->localtime() - (trigger_time + offset));
+  }
+
+  return StatusCode::SUCCESS;
 }
 
 // Find the baseline
@@ -248,185 +337,129 @@ WaveformReconstructionTool::findBaseline(const RawWaveform& raw_wave,
     // Save baseline to hit collection object
     hit->set_baseline_mean(raw_wave.mv_per_bit()*baseline.mean);
     hit->set_baseline_rms(raw_wave.mv_per_bit()*baseline.rms);
+    ATH_MSG_DEBUG("Baseline found with mean = " << hit->baseline_mean()
+		  << " mV and rms = " << hit->baseline_rms()
+		  << " mV");
   }
 
   return baseline;
 }
 
-StatusCode
-WaveformReconstructionTool::reconstruct(const RawWaveform& raw_wave,
-					const xAOD::WaveformClock* clock, 
-					xAOD::WaveformHitContainer* container) const {
-
-  ATH_MSG_DEBUG(" reconstruct called ");
-
-  // Check the container
-  if (!container) {
-    ATH_MSG_ERROR("WaveformHitCollection passed to reconstruct() is null!");
-    return StatusCode::FAILURE;
-  }
-
-  //
-  // We always want to create at least one hit, so create it here
-  xAOD::WaveformHit* hit = new xAOD::WaveformHit();
-  container->push_back(hit);
+// Fill the raw hit parameters
+void
+WaveformReconstructionTool::fillRawHitValues(const RawWaveform& wave,
+					   int lo_edge, int hi_edge, 
+					   xAOD::WaveformHit* hit) const {
 
-  // Set digitizer channel and identifier
-  hit->set_channel(raw_wave.channel());
-  hit->set_id(raw_wave.identify32().get_compact());
+  // First, make sure we don't overflow the waveform range
+  if (lo_edge < 0) lo_edge = 0;
+  if (hi_edge >= int(wave.size())) hi_edge = wave.size() - 1;
 
-  // 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() << "!");
+  ATH_MSG_DEBUG("Fill channel " << wave.channel() 
+		<< " waveform from sample " << lo_edge << " to " << hi_edge);
 
-    hit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_MISSING);
-    return StatusCode::SUCCESS;
-  } 
+  // Fill hit window with data from wave
+  std::vector<float> wtime(hi_edge-lo_edge+1);
+  std::vector<float> wwave(hi_edge-lo_edge+1);
 
-  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() 
-		     << " not equal to number of samples " << raw_wave.n_samples());
-
-    hit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_INVALID);
-    return StatusCode::SUCCESS;
+  for (int i=lo_edge; i<=hi_edge; i++) {
+    unsigned int j = i-lo_edge;
+    wtime[j] = 2.*i; // 2ns per sample at 500 MHz
+    wwave[j] = hit->baseline_mean() - wave.mv_per_bit() * wave.adc_counts()[i];
   }
 
-  // Find the baseline
-  WaveformBaselineData baseline = findBaseline(raw_wave, hit);
-
-  // 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
-  // Time in ns from start of readout
-  unsigned int size = raw_wave.adc_counts().size();
-  std::vector<float> time(size);  
-  for (unsigned int i=0; i<size; i++)
-    time[i] = 2.*i;
-
-  // Baseline subtracted (and inverted) ADC waveform values
-  std::vector<float> wave(raw_wave.adc_counts().begin(), raw_wave.adc_counts().end());
-  for (auto& element : wave)
-    element = baseline.mean - element;
-
-  bool first = true;
-
-  // Now we iteratively find peaks and fit
-  while(true) {
+  hit->set_time_vector(wtime);
+  hit->set_wave_vector(wwave);
 
-    //
-    // Find peak in array and return time and value arrays
-    // This range of data is also *removed* from original arrays
-    std::vector<float> wtime;
-    std::vector<float> wwave;
-
-    // All done if we don't have any peaks above threshold
-    // If we do find a significant peak, fill the window
-    if (! findPeak(baseline, time, wave, wtime, wwave) ) {
-      if (first) hit->set_status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED);
-      break;
-    }
-
-    //
-    // Create new hit to fill
-    if (!first) {
-      hit = new xAOD::WaveformHit();
-      container->push_back(hit);
-      hit->set_status_bit(xAOD::WaveformStatus::SECONDARY);
-    }
-    first = false;
-
-    //
-    // Save windowed waveform to Hit object
-    hit->set_channel(raw_wave.channel());
-    hit->set_baseline_mean(baseline.mean);
-    hit->set_baseline_rms(baseline.rms);
-    hit->set_time_vector(wtime);
-    hit->set_wave_vector(wwave);
+  // Set raw values
+  WaveformFitResult raw = findRawHitValues(wtime, wwave);
+  hit->set_peak(raw.peak);
+  hit->set_mean(raw.mean);
+  hit->set_width(raw.sigma);
+  hit->set_integral(raw.integral);
+  hit->set_localtime(raw.mean);
+  hit->set_raw_peak(raw.peak);
+  hit->set_raw_integral(raw.integral);
 
-    //
-    // Find some raw values
-    WaveformFitResult raw = findRawHitValues(wtime, wwave);
-    hit->set_peak(raw.peak);
-    hit->set_mean(raw.mean);
-    hit->set_width(raw.sigma);
-    hit->set_integral(raw.integral);
-    hit->set_localtime(raw.mean);
-    hit->set_raw_peak(raw.peak);
-    hit->set_raw_integral(raw.integral);
+}
 
-    //
-    // Perform Gaussian fit to waveform
-    WaveformFitResult gfit = fitGaussian(raw, wtime, wwave);
-    if (! gfit.valid) {
-      // Lets try again with a more restricted width
-      ATH_MSG_WARNING( " Gaussian waveform fit failed with width " << raw.sigma << " try reducing width to 1 " );
-      raw.sigma = 1.;
-      gfit = fitGaussian(raw, wtime, wwave);
-      if (!gfit.valid) {
-	hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED);
-      }
-    } 
+// Reconstruct a hit from the RawWaveform in the range specified
+// Range is in units digitizer samples (not ns)
+void
+WaveformReconstructionTool::reconstructHit(xAOD::WaveformHit* hit) const {
 
-    // Fit results (or raw if it failed)
-    hit->set_peak(gfit.peak);
-    hit->set_mean(gfit.mean);
-    hit->set_width(gfit.sigma);
-    hit->set_integral(gfit.integral);
-    hit->set_localtime(gfit.time);
+  // Time and waveform vectors
+  // Don't use reference as we may modify this below
+  std::vector<float> wtime = hit->time_vector();
+  std::vector<float> wwave = hit->wave_vector();
 
-    //
-    // Check for overflow
-    if (m_removeOverflow && findOverflow(baseline, wtime, wwave)) {
-      ATH_MSG_INFO("Found waveform overflow");
-      hit->set_status_bit(xAOD::WaveformStatus::WAVE_OVERFLOW);
-    }
+  ATH_MSG_DEBUG("Reconstruct channel " << hit->channel() 
+		<< " waveform from " << wtime.front() 
+		<< " to " << wtime.back());
 
-    //
-    // Perform CB fit
-    WaveformFitResult cbfit = fitCBall(gfit, wtime, wwave);
-    if (! cbfit.valid) {
-      ATH_MSG_WARNING("CrystalBall fit failed!");
-      // Still have gaussian parameters as an estimate
-      hit->set_status_bit(xAOD::WaveformStatus::CBFIT_FAILED);
-    } else {
-      hit->set_peak(cbfit.peak);
-      hit->set_mean(cbfit.mean);
-      hit->set_width(cbfit.sigma);
-      hit->set_integral(cbfit.integral);
-      hit->set_localtime(cbfit.time);
-
-      hit->set_alpha(cbfit.alpha);
-      hit->set_nval(cbfit.nval);
-    }
+  // Fill values needed for fit (peak, mean, and sigma)
+  WaveformFitResult raw;
+  raw.peak = hit->peak();
+  raw.mean = hit->mean();
+  raw.sigma = hit->width();
 
-    //
-    // Find time from clock
-    if (!clock || (clock->frequency() <= 0.)) {
-      hit->set_status_bit(xAOD::WaveformStatus::CLOCK_INVALID);
-      hit->set_bcid_time(-1.);
-    } else {
-      hit->set_bcid_time(clock->time_from_clock(hit->localtime()));
+  //
+  // Perform Gaussian fit to waveform
+  WaveformFitResult gfit = fitGaussian(raw, wtime, wwave);
+  if (! gfit.valid) {
+    // Lets try again with a more restricted width
+    ATH_MSG_WARNING( " Gaussian waveform fit failed with width " << raw.sigma << " try reducing width to 1 " );
+    raw.sigma = 1.;
+    gfit = fitGaussian(raw, wtime, wwave);
+    if (!gfit.valid) {
+      hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED);
     }
+  } 
 
-    if (! m_findMultipleHits) break;
-
-  } // End of loop over waveform data
+  // Fit results (or raw if it failed)
+  hit->set_peak(gfit.peak);
+  hit->set_mean(gfit.mean);
+  hit->set_width(gfit.sigma);
+  hit->set_integral(gfit.integral);
+  hit->set_localtime(gfit.time);
+  
+  //
+  // Check for overflow
+  if (m_removeOverflow && findOverflow(hit->baseline_mean(), wtime, wwave)) {
+    ATH_MSG_INFO("Found waveform overflow");
+    hit->set_status_bit(xAOD::WaveformStatus::WAVE_OVERFLOW);
+  }
+  
+  //
+  // Perform CB fit
+  WaveformFitResult cbfit = fitCBall(gfit, wtime, wwave);
+  if (! cbfit.valid) {
+    ATH_MSG_WARNING("CrystalBall fit failed for channel " << hit->channel() << "!");
+    // Still have gaussian parameters as an estimate
+    hit->set_status_bit(xAOD::WaveformStatus::CBFIT_FAILED);
+  } else {
+    hit->set_peak(cbfit.peak);
+    hit->set_mean(cbfit.mean);
+    hit->set_width(cbfit.sigma);
+    hit->set_integral(cbfit.integral);
+    hit->set_localtime(cbfit.time);
+    
+    hit->set_alpha(cbfit.alpha);
+    hit->set_nval(cbfit.nval);
+  }
 
-  ATH_MSG_DEBUG( "WaveformReconstructionTool finished for channel " 
-		 << raw_wave.channel() << " container size= " << container->size());
+  ATH_MSG_DEBUG("Done reconstructing channel " << hit->channel() 
+		<< " waveform from " << wtime.front() << " to " << wtime.back());
 
-  return StatusCode::SUCCESS;
 }
 
-bool
+// Returns location of peak in array wave
+// Return value is -1 if peak is below threshold
+int
 WaveformReconstructionTool::findPeak(WaveformBaselineData& baseline, 
-   std::vector<float>& time, std::vector<float>& wave,
-   std::vector<float>& windowed_time, std::vector<float>& windowed_wave) const {
+				     float threshold,
+				     std::vector<float>& wave) const 
+{
 
   ATH_MSG_DEBUG("findPeak called");
 
@@ -436,44 +469,31 @@ WaveformReconstructionTool::findPeak(WaveformBaselineData& baseline,
   ATH_MSG_DEBUG( "Found peak value " << maxval << " at position " << imax );
 
   // Check if this is over threshold (in sigma)
-  if (maxval < m_peakThreshold*baseline.rms) {
+  if (maxval < threshold*baseline.rms) {
     ATH_MSG_DEBUG("Failed threshold");
-    return false;
+    return -1;
   }
 
-  // Make a window around this peak, values are in bins, so units of 2ns
-  // Ensure our window is within the vector range
-  int lo_edge = ((int(imax) + m_windowStart) >= 0 ? (imax + m_windowStart) : 0); 
-  int hi_edge = ((imax + m_windowStart + m_windowWidth) < wave.size() ? (imax + m_windowStart + m_windowWidth) : wave.size());
-  
-  ATH_MSG_DEBUG("Windowing waveform from " << lo_edge << " to " << hi_edge);
-  windowed_time = std::vector<float> (time.begin()+lo_edge, time.begin()+hi_edge);
-  windowed_wave = std::vector<float> (wave.begin()+lo_edge, wave.begin()+hi_edge);
-
-  // Remove these values from the original arrays so we can iterate
-  time.erase(time.begin()+lo_edge, time.begin()+hi_edge);
-  wave.erase(wave.begin()+lo_edge, wave.begin()+hi_edge);
-
-  return true;
+  return imax;
 }
 
 bool
-WaveformReconstructionTool::findOverflow(const WaveformBaselineData& base, 
+WaveformReconstructionTool::findOverflow(float baseline,
 	      std::vector<float>& time, std::vector<float>& wave) const {
 
   auto peakloc = std::max_element(wave.begin(), wave.end());
 
   // If peak value is less than baseline, we have no overflow
-  if (*peakloc < int(base.mean)) return false;
+  if (*peakloc < baseline) return false;
 
   ATH_MSG_DEBUG("Removing overflows from waveform with length " << wave.size());
 
   // We have an overflow, remove all elements that are overflowing
   unsigned int i = peakloc - wave.begin();
   for (; i<wave.size(); i++) {
-    if (wave[i] < int(base.mean)) continue;
+    if (wave[i] < baseline) continue;
 
-    ATH_MSG_DEBUG("Removing position "<< i<< " with value " << wave[i] << " > " << int(base.mean));
+    ATH_MSG_DEBUG("Removing position "<< i<< " with value " << wave[i] << " > " << baseline);
     // This is an overflow, remove elements
     time.erase(time.begin() + i);
     wave.erase(wave.begin() + i);
@@ -772,7 +792,7 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit,
   TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", "");
 
   if (!cbfitptr->IsValid()) {
-    ATH_MSG_WARNING( " First Crystal Ball waveform fit failed! ");
+    ATH_MSG_DEBUG( " First Crystal Ball waveform fit failed! ");
   }
 
   // Now try releasing the tail parameter
@@ -787,7 +807,7 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit,
   cbfit.valid = (cbfit.fit_status == 0);
 
   if (!cbfitptr->IsValid()) {
-    ATH_MSG_WARNING( " Crystal Ball waveform fit failed! ");
+    ATH_MSG_DEBUG( " Full Crystal Ball waveform fit failed! ");
   } else {
     // Improve estimation with fit results
     cbfit.peak = cbfitptr->Parameter(0);
diff --git a/Waveform/WaveRecTools/src/WaveformReconstructionTool.h b/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
index 7a39f883434051173099d3702e9643fb64cc9f51..b23b60e9491ff2989601c436f4723ca7fa76dfc2 100644
--- a/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
+++ b/Waveform/WaveRecTools/src/WaveformReconstructionTool.h
@@ -20,6 +20,9 @@
 #include "WaveformBaselineData.h"
 #include "WaveformFitResult.h"
 
+// Tool classes
+#include "WaveformConditionsTools/IWaveformTimingTool.h"
+
 //Gaudi
 #include "GaudiKernel/ToolHandle.h"
 
@@ -37,16 +40,18 @@ 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 primary hits from waveform (in trigger window)
+  virtual StatusCode reconstructPrimary(const RawWaveform& wave,
+					xAOD::WaveformHitContainer* hitContainer) const;
+
+  /// Reconstruct primary hits from waveform (in trigger window)
+  virtual StatusCode reconstructSecondary(const RawWaveform& wave,
+					xAOD::WaveformHitContainer* hitContainer) const;
+
+  /// Set local hit times from LHC clock
+  virtual StatusCode setLocalTime(const xAOD::WaveformClock* clock,
+				  xAOD::WaveformHitContainer* container) const;
 
-  /// Reconstruct hits from waveform
-  
-  virtual StatusCode reconstruct(const RawWaveform& wave,
-				 const xAOD::WaveformClock* clock,
-				 xAOD::WaveformHitContainer* hitContainer) const;
 
  private:
 
@@ -54,6 +59,9 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
   // Baseline Estimation Parameters
   BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", false};
 
+  ToolHandle<IWaveformTimingTool> m_timingTool
+    {this, "WaveformTimingTool", "WaveformTimingTool"};
+
   // Minimum number of samples needed to calculate simple baseline
   // Just average these first n values
   IntegerProperty m_samplesForBaselineAverage{this, "SamplesForBaselineAverage", 40};
@@ -79,32 +87,40 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
   FloatProperty m_baselineFitWindow{this, "BaselineFitWindow", 2.};
 
   //
-  // Peak threshold (in sigma of baseline RMS) to find a hit
-  FloatProperty m_peakThreshold{this, "PeakThreshold", 10.};
+  // Peak threshold (in sigma of baseline RMS)
+  // Primary threshold is requirement to try a fit for the in-time window
+  // Secondary threshold is requirement to produce a secondary hit
+  // from a local maximum
+  FloatProperty m_primaryPeakThreshold{this, "PrimaryPeakThreshold", 5.};
+  FloatProperty m_secondaryPeakThreshold{this, "SecondaryPeakThreshold", 10.};
  
   //
   // Window to define fitting range, in samples (2ns/sample)
-  IntegerProperty m_windowStart{this, "FitWindowStart", -15};
+  IntegerProperty m_windowStart{this, "FitWindowStart", -20};  
   IntegerProperty m_windowWidth{this, "FitWindowWidth", 60};
 
   //
   // Remove overflow values from CB fit
   BooleanProperty m_removeOverflow{this, "RemoveOverflow", true};
 
-  //
-  // Look for more than one hit in each channel
-  BooleanProperty m_findMultipleHits{this, "FindMultipleHits", false};
-
   //
   // Fraction of peak to set local hit time
-  FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.45};
+  FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.5};
 
   //
-  // 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.};
+  // When looking for secondary hits with a primary found above threshold
+  // should we look before or after the primary hit?
+  BooleanProperty m_findSecondaryBefore{this, "FindSecondaryBefore", true};
+  BooleanProperty m_findSecondaryAfter{this, "FindSecondaryAfter", false};
+
+  // Reco algorithms
+  // Fill hit with raw data from waveform
+  void fillRawHitValues(const RawWaveform& wave,
+			int lo_edge, int hi_edge,
+			xAOD::WaveformHit* hit) const;
+
+  // Perform fits to WaveformHit data
+  void reconstructHit(xAOD::WaveformHit* hit) const;
 
   // Baseline algorithms
   WaveformBaselineData& findSimpleBaseline(const RawWaveform& wave) const;
@@ -113,12 +129,10 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
 				     xAOD::WaveformHit* hit) const;
 
 
-  // Find peak in wave, return windowed region in windowed_time and windowed_wave
-  // Windowed region is removed from original vectors
-  // Returns true if peak found, false if not
-  bool findPeak(WaveformBaselineData& baseline,
-		std::vector<float>& time, std::vector<float>& wave,
-		std::vector<float>& windowed_time, std::vector<float>& windowed_wave) const;
+  // Find peak in wave, return index to peak position, or -1 if 
+  // peak isn't greater than threshold
+  int findPeak(WaveformBaselineData& baseline, float threshold, 
+		std::vector<float>& wave) const;
 
   // Get estimate from waveform data itself
   WaveformFitResult& findRawHitValues(const std::vector<float> time, 
@@ -130,7 +144,7 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
 				 const std::vector<float> wave) const;
 
   // Find overflows and remove points from arrays
-  bool findOverflow(const WaveformBaselineData& baseline, 
+  bool findOverflow(float baseline, 
 		    std::vector<float>& time, std::vector<float>& wave) const;
 
   // Fit windowed data to CrystalBall function
@@ -138,12 +152,6 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
 			      const std::vector<float> time, 
 			      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/Waveform/WaveformConditions/WaveCondUtils/CMakeLists.txt b/Waveform/WaveformConditions/WaveCondUtils/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1c295bd9cab84b213f7ec0227a3cb4f95486edc8
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveCondUtils/CMakeLists.txt
@@ -0,0 +1,9 @@
+################################################################################
+# Package: WaveCondUtils
+################################################################################
+
+# Declare the package name:
+atlas_subdir( WaveCondUtils )
+
+atlas_install_scripts( scripts/*.sh scripts/*.py )
+
diff --git a/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py
new file mode 100755
index 0000000000000000000000000000000000000000..b3a1c64a23ad23ae70eaa000e3cbd76064bebfcd
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveCondUtils/scripts/makeTimingDB.py
@@ -0,0 +1,221 @@
+#!/bin/env python
+
+# Requires python 3.8 or higher
+#
+# Can test results with
+# AtlCoolConsole.py "sqlite://;schema=waveform_reco.db;dbname=OFLP200"
+
+filename = 'waveform_reco.db'
+
+# Nominal trigger time in ns
+nominal_data = {
+    0: 820.,
+    4272: 830.,
+    6525: 820.
+}
+
+offset_channels = 16
+
+# Run 
+# 0 - initial data
+# 3395 - Testbeam
+# 
+
+ehn1_offsets = [ -20., -20., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ]
+ti12_offsets = [ -20., -20., -20., -20., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ]
+
+offset_data = {
+  0:    [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ],
+# Initial TI12
+  1324: [ -10., -10., -10., -10.,   0.,   0., 0., 0., 0., 0., 0., 0., 18., 18., 0., 0. ],
+# Testbeam geometry
+  3247: [ -10., -10., -10., -10., -10., -10., 15., 15., -20., -20., 0., 0., 0., 0., 0., 0. ],
+# TI12
+  4272: ti12_offsets,
+# EHN1 (interleaved with TI12 running)
+  4360: ehn1_offsets,
+  4399: ti12_offsets,
+  4409: ehn1_offsets,
+  4411: ti12_offsets,
+  4429: ehn1_offsets,
+  4439: ti12_offsets,
+  4876: ehn1_offsets,
+  4892: ti12_offsets,
+  4904: ehn1_offsets,
+  4912: ti12_offsets,
+  4954: ehn1_offsets,
+  4989: ti12_offsets,
+  4991: ehn1_offsets,
+  4993: ti12_offsets,
+  4996: ehn1_offsets,
+  4997: ti12_offsets,
+  5042: ehn1_offsets,
+  5050: ti12_offsets,
+# IFT and VetoNu installed
+  6525: [ -10., -10., -10., -10., -25., -25., 0., 0., 0., 0., 0., 0., 18., 18., 0., 0. ]
+}
+
+attr_list_desc = '<timeStamp>run-lumi</timeStamp><addrHeader><address_header service_type="71" clid="40774348" /></addrHeader><typeName>AthenaAttributeList</typeName>'
+
+cond_attr_list_desc = '<timeStamp>run-lumi</timeStamp><addrHeader><address_header clid="1238547719" service_type="71" /></addrHeader><typeName>CondAttrListCollection</typeName>'
+
+maxInt32 = 0xFFFFFFFF
+
+
+# Look for data entry errors
+
+print('Validating nominal data')
+
+lastRun = -1
+for run, data in nominal_data.items():
+    assert isinstance(run, int), 'Run number is not integer'
+    assert isinstance(data, float), 'Time is not float'
+    assert run > lastRun, 'Run numbers out of order'
+    assert run <= maxInt32, 'Run number out of range'
+    lastRun = run
+
+print('Validating offset data')
+lastRun = -1
+for run, data in offset_data.items():
+    assert isinstance(run, int), 'Run number is not integer'
+    assert run > lastRun, 'Run numbers out of order'
+    assert run <= maxInt32, 'Run number out of range'
+    lastRun = run
+    assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
+    for i in range(offset_channels):
+        assert isinstance(data[i], float), 'Offset time is not float'
+
+# Data looks OK
+
+
+from PyCool import cool
+
+dbSvc = cool.DatabaseSvcFactory.databaseService()
+connectString = f'sqlite://;schema={filename};dbname=CONDBR3'
+
+print('Creating database')
+
+dbSvc.dropDatabase( connectString )
+db = dbSvc.createDatabase( connectString )
+
+# Nominal trigger times
+nominalSpec = cool.RecordSpecification()
+nominalSpec.extend( 'NominalTriggerTime', cool.StorageType.Float )
+
+nominalFolderSpec = cool.FolderSpecification(cool.FolderVersioning.SINGLE_VERSION, nominalSpec)
+nominalFolder = db.createFolder('/WAVE/DAQ/Timing', nominalFolderSpec, attr_list_desc, True)
+
+# There should be one record entered per IOV
+lastValid = cool.ValidityKeyMax
+for firstValidRun, time in reversed(nominal_data.items()):
+    firstValid = (firstValidRun << 32)
+    nominalRecord = cool.Record(nominalSpec)
+    nominalRecord[ 'NominalTriggerTime' ] = float(time)
+    nominalFolder.storeObject( firstValid, lastValid, nominalRecord, cool.ChannelId(0))
+    lastValid = ((firstValidRun - 1) << 32) | (cool.ValidityKeyMax & 0x00000000FFFFFFFF)
+
+
+# Trigger offset times
+
+offsetSpec = cool.RecordSpecification()
+offsetSpec.extend( 'TriggerOffset', cool.StorageType.Float )
+
+offsetFolderSpec = cool.FolderSpecification(cool.FolderVersioning.SINGLE_VERSION, offsetSpec)
+offsetFolder = db.createFolder('/WAVE/DAQ/TimingOffset', offsetFolderSpec, cond_attr_list_desc, True)
+
+# There should be one record entered per IOV
+lastValid = cool.ValidityKeyMax
+for firstValidRun, offset_list in reversed(offset_data.items()):
+    firstValid = (firstValidRun << 32)
+    for channel in range(offset_channels):
+        offsetRecord = cool.Record(offsetSpec)
+        offsetRecord[ 'TriggerOffset' ] = float(offset_list[channel])
+        offsetFolder.storeObject( firstValid, lastValid, offsetRecord, cool.ChannelId(channel) )
+
+    lastValid = ((firstValidRun - 1) << 32) | (cool.ValidityKeyMax & 0x00000000FFFFFFFF)
+
+
+db.closeDatabase()
+
+print('Database completed')
+
+print('Working on MC database')
+
+# Nominal data 
+nominal_data = {
+    0: 820.
+}
+# No offsets by default
+offset_data = {
+  0: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
+}
+
+# Validate again just in case
+print('Validating nominal data')
+
+lastRun = -1
+for run, data in nominal_data.items():
+    assert isinstance(run, int), 'Run number is not integer'
+    assert isinstance(data, float), 'Time is not float'
+    assert run > lastRun, 'Run numbers out of order'
+    assert run <= maxInt32, 'Run number out of range'
+    lastRun = run
+
+print('Validating offset data')
+lastRun = -1
+for run, data in offset_data.items():
+    assert isinstance(run, int), 'Run number is not integer'
+    assert run > lastRun, 'Run numbers out of order'
+    assert run <= maxInt32, 'Run number out of range'
+    lastRun = run
+    assert len(data) == offset_channels, 'Offset data does not have '+str(offset_channels)+' entries'
+    for i in range(offset_channels):
+        assert isinstance(data[i], float), 'Offset time is not float'
+
+# Data looks OK
+
+connectString = f'sqlite://;schema={filename};dbname=OFLP200'
+
+dbSvc.dropDatabase( connectString )
+db = dbSvc.createDatabase( connectString )
+
+# Nominal trigger times
+nominalSpec = cool.RecordSpecification()
+nominalSpec.extend( 'NominalTriggerTime', cool.StorageType.Float )
+
+nominalFolderSpec = cool.FolderSpecification(cool.FolderVersioning.SINGLE_VERSION, nominalSpec)
+nominalFolder = db.createFolder('/WAVE/DAQ/Timing', nominalFolderSpec, attr_list_desc, True)
+
+# There should be one record entered per IOV
+lastValid = cool.ValidityKeyMax
+for firstValidRun, time in reversed(nominal_data.items()):
+    firstValid = (firstValidRun << 32)
+    nominalRecord = cool.Record(nominalSpec)
+    nominalRecord[ 'NominalTriggerTime' ] = float(time)
+    nominalFolder.storeObject( firstValid, lastValid, nominalRecord, cool.ChannelId(0))
+    lastValid = ((firstValidRun - 1) << 32) | (cool.ValidityKeyMax & 0x00000000FFFFFFFF)
+
+
+# Trigger offset times
+
+offsetSpec = cool.RecordSpecification()
+offsetSpec.extend( 'TriggerOffset', cool.StorageType.Float )
+
+offsetFolderSpec = cool.FolderSpecification(cool.FolderVersioning.SINGLE_VERSION, offsetSpec)
+offsetFolder = db.createFolder('/WAVE/DAQ/TimingOffset', offsetFolderSpec, cond_attr_list_desc, True)
+
+# There should be one record entered per IOV
+lastValid = cool.ValidityKeyMax
+for firstValidRun, offset_list in reversed(offset_data.items()):
+    firstValid = (firstValidRun << 32)
+    for channel in range(offset_channels):
+        offsetRecord = cool.Record(offsetSpec)
+        offsetRecord[ 'TriggerOffset' ] = float(offset_list[channel])
+        offsetFolder.storeObject( firstValid, lastValid, offsetRecord, cool.ChannelId(channel) )
+
+    lastValid = ((firstValidRun - 1) << 32) | (cool.ValidityKeyMax & 0x00000000FFFFFFFF)
+
+
+db.closeDatabase()
+
+print('Database completed')
diff --git a/Waveform/WaveformConditions/WaveCondUtils/scripts/wave_timing_check.py b/Waveform/WaveformConditions/WaveCondUtils/scripts/wave_timing_check.py
new file mode 100755
index 0000000000000000000000000000000000000000..6c8f6332ff42e1ffc06fb168d290539c1db89208
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveCondUtils/scripts/wave_timing_check.py
@@ -0,0 +1,337 @@
+#!/usr/bin/env python3
+#
+import os
+import sys
+import math
+import array
+import itertools
+
+# Triggers: 0x01 - calo, 0x02 - veto, 0x03 - timing, 0x10 - random
+
+def usage():
+    print("Usage: timing_check.py <filename>|<dirname> [triggermask]")
+
+if len(sys.argv) == 1:
+    usage()
+    sys.exit(-1)
+
+# Extract tracker station requirements
+if len(sys.argv) == 3:
+    trigmask = int(sys.argv[2])
+    extra = f"_{triggermask}"
+else:
+    trigmask = 0xFF
+    extra = ''
+
+from pathlib import Path
+
+import ROOT
+ROOT.xAOD.Init().ignore()
+ROOT.xAOD.AuxContainerBase()
+os.environ["XAOD_ACCESSTRACER_FRACTION"] = "0.0"
+
+#
+# Open file or files
+pathname = Path(sys.argv[1])
+
+# Is this a directory?
+if pathname.is_dir():
+    print(f"Opening files in directory {pathname.name}")
+
+    t2 = ROOT.TChain("CollectionTree")
+    nfiles = t2.Add(str(pathname)+'/Faser-Physics*.root')
+
+    if (nfiles == 0):
+        print(f"TChain found no files!")
+        usage()
+        sys.exit(0)
+
+    # Make transient tree
+    t1 = ROOT.xAOD.MakeTransientTree(t2)
+
+    # Make output file name
+    outfile = pathname.name + "_timing"+extra+".pdf"
+
+    print(f"TChain found {nfiles} files with {t2.GetEntries()} events")
+
+    avperfile = t2.GetEntries() / nfiles
+
+# Is this a file?
+elif pathname.is_file():
+    print(f"Opening file {pathname.name}")
+
+    t2 = ROOT.TChain("CollectionTree")
+    nfiles = t2.Add(str(pathname))
+
+    if (nfiles != 1):
+        print(f"TChain error opening file!")
+        usage()
+        sys.exit(0)
+
+    print(f"Opened file with {t2.GetEntries()} events")
+
+    avperfile = t2.GetEntries()
+
+    # Make transient tree
+    t1 = ROOT.xAOD.MakeTransientTree(t2)
+
+    # Make outfile name from input
+    outfile = pathname.stem + "_timing"+extra+".pdf"
+
+# Neither?
+else:
+    print(f"Can't understand {pathname.name}")
+    usage()
+    sys.exit(-1)
+
+class ClockPlots:
+
+    def __init__(self):
+
+        # Ranges for plots
+        self.freq_bins = 80
+        self.freq_lo = 40.0
+        self.freq_hi = 40.2
+
+        self.th_bins = 100
+
+    def init(self, tree):
+
+        self.h_freq = ROOT.TH1I("", "Clock Frequency", self.freq_bins, self.freq_lo, self.freq_hi)
+        self.h_freq.GetXaxis().SetTitle("Clock Frequency (MHz)")
+        self.h_freq.GetYaxis().SetTitle("Events")
+        #self.h_freq.Sumw2()
+
+        self.h_phase = ROOT.TH1I("", "Clock Phase", 60, 2*(-3.1416), 2*3.1416)
+        self.h_phase.GetXaxis().SetTitle("Clock Phase")
+        self.h_phase.GetYaxis().SetTitle("Events")
+
+        self.h_amp = ROOT.TH1I("", "Amplitude", 50, 0, 2000.)
+        self.h_amp.GetXaxis().SetTitle("Clock Amplitude (mV)")
+        self.h_amp.GetYaxis().SetTitle("Events")
+
+        self.h_off = ROOT.TH1I("", "Offset", 50, 0, 2000.)
+        self.h_off.GetXaxis().SetTitle("Clock Offset (mV)")
+        self.h_off.GetYaxis().SetTitle("Events")
+
+    def fill(self, tree):
+
+        # First, create the histograms
+        self.init(tree)
+
+        # Iterate over all entries
+        nev = tree.GetEntries()
+        iev = 0
+        for ev in tree:
+            self.h_freq.Fill(ev.WaveformClock.frequency())
+            self.h_phase.Fill(ev.WaveformClock.phase())
+            self.h_amp.Fill(ev.WaveformClock.amplitude())
+            self.h_off.Fill(ev.WaveformClock.dc_offset())
+
+            # Protect against reading off the end
+            iev += 1
+            if iev == nev: break
+
+    def draw(self, canvas, outfile):
+
+        # Under/overflows, mean, rms, and entries
+        ROOT.gStyle.SetOptStat(111110)
+
+        canvas.Clear()
+        canvas.Divide(2,2)
+        canvas.cd(1)
+        self.h_freq.Draw()
+        canvas.cd(2)
+        self.h_phase.Draw()
+        canvas.cd(3)
+        self.h_amp.Draw()
+        canvas.cd(4)
+        self.h_off.Draw()
+        canvas.Update()
+        canvas.Print(outfile)
+
+    def print_stats(self):
+
+        freq_mean = self.h_freq.GetMean()
+        freq_rms = self.h_freq.GetStdDev()
+        freq_n = self.h_freq.GetEntries()
+        print(f"LHC Clock: {freq_mean:.6} +/- {freq_rms/math.sqrt(freq_n):.6}")
+
+class WavePlots:
+
+    def __init__(self, triggerMask=0xFF):
+
+        # Number of waveforms channels
+        self.nchan = 15
+
+        # Trigger mask
+        self.mask = triggerMask
+
+        self.chan_hist_list = []
+        self.log_list = []
+
+        # Maaximum peak value
+        self.peak_max = 16000.
+
+    def init(self, tree):
+
+        # Keyed by channel
+        self.createChannelHist('h_localtime', 40, 750, 950, "Local Time")
+        self.createChannelHist('h_triggertime', 40, -80, 80, "Trigger Time")
+        self.createChannelHist('h_bcidtime', 50, -10, 40, "BCID Time")
+
+    def createChannelHist(self, name, nbins, xlo, xhi, xtitle='', ytitle='Waveforms', stats=True, log=False):
+
+        setattr(self, name, dict())
+        x = getattr(self, name)
+        for chan in range(self.nchan):
+            x[chan] = ROOT.TH1I("", "", nbins, xlo, xhi)
+            if len(xtitle) > 0:
+                x[chan].GetXaxis().SetTitle(f"Ch {chan} {xtitle}")
+            if len(ytitle) > 0:
+                x[chan].GetYaxis().SetTitle(ytitle)
+            x[chan].SetStats(stats)
+
+        self.chan_hist_list.append(name)
+        if log:
+            self.log_list.append(name)
+
+    def fill(self, tree):
+
+        # First, create the histograms
+        self.init(tree)
+
+        # Iterate over all entries
+        nev = tree.GetEntries()
+        iev = 0
+        for ev in tree:
+
+            time = ev.EventInfo.timeStamp()
+            trig = ev.FaserTriggerData.tap()
+
+            if not (trig & self.mask): 
+                iev += 1
+                if iev == nev: 
+                    break
+                else: 
+                    continue
+
+            # Process waveforms
+            try:
+                wave_list = itertools.chain(ev.CaloWaveformHits, ev.PreshowerWaveformHits, ev.TriggerWaveformHits, ev.VetoWaveformHits, ev.VetoNuWaveformHits)
+            except:
+                wave_list = itertools.chain(ev.CaloWaveformHits, ev.PreshowerWaveformHits)
+
+            for wave in wave_list:
+
+                channel = wave.channel()
+
+                # Check if failed threshold
+                if wave.status_bit(0): continue
+
+                # Fill fit parameters
+                self.h_localtime[channel].Fill(wave.localtime())
+                self.h_triggertime[channel].Fill(wave.trigger_time())
+                self.h_bcidtime[channel].Fill(wave.bcid_time())
+
+            # End of loop over waveforms
+
+            # Protect against reading off the end
+            iev+=1
+            if iev == nev: break
+
+        # End of loop over events
+
+        # Put overflows in last bin of plots
+        self.fixOverflow(self.h_localtime)
+        self.fixOverflow(self.h_triggertime)
+
+    def fixOverflow(self, hdict):
+
+        for h in hdict.values():
+
+            if h.GetNbinsY() == 1:
+                self.fixOverflow1D(h)
+            else:
+                self.fixOverflow2D(h)
+
+    def fixOverflow1D(self, hist):
+        nbins = hist.GetNbinsX()
+        nlast = hist.GetBinContent(nbins)
+        nover = hist.GetBinContent(nbins+1)
+        hist.SetBinContent(nbins, nlast+nover)
+
+    def fixOverflow2D(self, hist):
+        nbx = hist.GetNbinsX()
+        nby = hist.GetNbinsY()
+
+        for ibinx in range(nbx+1):
+            nlast = hist.GetBinContent(ibinx, nby)
+            nover = hist.GetBinContent(ibinx, nby+1)
+            hist.SetBinContent(ibinx, nby, nlast+nover)
+
+        for ibiny in range(nby+1):
+            nlast = hist.GetBinContent(nbx, ibiny)
+            nover = hist.GetBinContent(nbx+1, ibiny)
+            hist.SetBinContent(nbx, ibiny, nlast+nover)
+
+        # Also the double overflow
+        nlast = hist.GetBinContent(nbx, nby)
+        nover = hist.GetBinContent(nbx+1, nby+1)
+        hist.SetBinContent(nbx, nby, nlast+nover)
+
+
+    def draw(self, canvas, outfile):
+
+        #
+        # Plot channel plots
+        for name in self.chan_hist_list:
+            canvas.Clear()
+            canvas.Divide(4,4)
+
+            if name in self.log_list:
+                setlog = True
+            else:
+                setlog = False
+
+            for chan in range(self.nchan):
+                canvas.cd(chan+1)
+                x = getattr(self, name)
+                x[chan].Draw()
+                if setlog: 
+                    ROOT.gPad.SetLogy(True)
+                else:
+                    ROOT.gPad.SetLogy(False)
+
+            canvas.Print(outfile)
+
+    def print_stats(self):
+
+        for chan in range(self.nchan):
+            local_mean = self.h_localtime[chan].GetMean()
+            trig_mean  = self.h_triggertime[chan].GetMean()
+            bcid_mean  = self.h_bcidtime[chan].GetMean()
+            print(f"Chan {chan:2}: Entries {int(self.h_localtime[chan].GetEntries()):8} Local {local_mean:6.1f} Trigger {trig_mean:6.2f} BCID {bcid_mean:6.2f}")
+
+#print("xAOD tree")
+#t1.Print()
+#print("non xAOD tree")
+#t2.Print()
+
+cp = ClockPlots()
+cp.fill(t1)
+
+# Triggers: 0x01 - calo, 0x02 - veto, 0x03 - timing, 0x10 - random
+wp = WavePlots(triggerMask=trigmask)
+wp.fill(t1)
+
+c = ROOT.TCanvas()
+c.Print(outfile+"[")
+
+cp.draw(c, outfile)
+wp.draw(c, outfile)
+
+c.Print(outfile+"]")
+
+cp.print_stats()
+wp.print_stats()
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformCableMappingTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformCableMappingTool.h
index 596d2bd8ce45189063f4952f7630ca0de0847c3e..461c087e0fb3957e90c1a1e50723e87a2baff996 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformCableMappingTool.h
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformCableMappingTool.h
@@ -38,6 +38,10 @@ class IWaveformCableMappingTool: virtual public IAlgTool {
   virtual WaveformCableMap getCableMapping(const EventContext& ctx) const = 0;
   virtual WaveformCableMap getCableMapping(void) const = 0;
 
+  virtual int getChannelMapping(const EventContext& ctx, const Identifier id) const = 0;
+  virtual int getChannelMapping(const Identifier id) const = 0;
+
+
 };
 
 //---------------------------------------------------------------------- 
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformTimingTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformTimingTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..de2a2dcbd1dde01c8a48e2a6a67822efc13e7519
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformTimingTool.h
@@ -0,0 +1,50 @@
+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS and FAsER collaborations
+*/
+
+/** @file IWaveformTimingTool.h Interface file for WaveformTimingTool.
+ * 
+ * Provides times and offsets (in ns) for different channels in the
+ * waveform digitizer.  This aligns the input signals for different
+ * path lengths and cable delays.
+ *
+ */
+
+// Multiple inclusion protection
+#ifndef IWAVEFORMTIMINGTOOL
+#define IWAVEFORMTIMINGTOOL
+
+//STL includes
+#include <map>
+
+//Gaudi Includes
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/EventContext.h"
+
+
+class IWaveformTimingTool: virtual public IAlgTool {
+
+ public:
+  
+  //----------Public Member Functions----------//
+  // Structors
+  virtual ~IWaveformTimingTool() = default; //!< Destructor
+
+  /// Creates the InterfaceID and interfaceID() method
+  DeclareInterfaceID(IWaveformTimingTool, 1, 0);
+
+  // Methods to return timing data
+
+  // Nominal trigger time (in ns) in the digitizer readout
+  virtual float nominalTriggerTime(void) const = 0;
+  virtual float nominalTriggerTime(const EventContext& ctx) const = 0;
+
+  // Channel-by-channel corrections to the nominal trigger time (in ns)
+  // A given channel should be centered at nominal + offset
+  virtual float triggerTimeOffset(int channel) const = 0;
+  virtual float triggerTimeOffset(const EventContext& ctx, int channel) const = 0;
+
+};
+
+//---------------------------------------------------------------------- 
+#endif // WAVEFORMTIMINGTOOL
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/python/WaveformTimingConfig.py b/Waveform/WaveformConditions/WaveformConditionsTools/python/WaveformTimingConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b95ed388208d7d4ab51a3d5621ab868d438c433
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/python/WaveformTimingConfig.py
@@ -0,0 +1,32 @@
+""" Define methods to configure WaveformTimingTool
+
+Copyright (C) 2022 CERN for the benefit of the FASER collaboration
+"""
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from IOVDbSvc.IOVDbSvcConfig import addFolders
+WaveformTimingTool=CompFactory.WaveformTimingTool
+
+def WaveformTimingToolCfg(flags, name="WaveformTimingTool", **kwargs):
+    """ Return a configured WaveformTimingTool"""
+    return WaveformTimingTool(name, **kwargs)
+
+def WaveformTimingCfg(flags, **kwargs):
+    """ Return configured ComponentAccumulator and tool for Waveform Timing 
+
+    WaveformTimingTool may be provided in kwargs
+    """
+
+    acc = ComponentAccumulator()
+    # tool = kwargs.get("WaveformTimingTool", WaveformTimingTool(flags))
+    # Probably need to figure this out!
+    dbInstance = kwargs.get("dbInstance", "TRIGGER_OFL")
+    if flags.Input.isMC:
+        dbname = "OFLP200"
+    else:
+        dbname = "CONDBR3"
+
+    acc.merge(addFolders(flags, "/WAVE/DAQ/Timing", dbInstance, className="AthenaAttributeList", db=dbname))
+    acc.merge(addFolders(flags, "/WAVE/DAQ/TimingOffset", dbInstance, className="CondAttrListCollection", db=dbname))
+    return acc
+
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
index 102aae0e79bda0b19b22b1969dad0103bf748d0f..530857354d1e1e82f5c522a72e8b6e0c7d6b5269 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.cxx
@@ -139,5 +139,98 @@ WaveformCableMappingTool::getCableMapping(void) const {
   return getCableMapping(ctx);
 }
 
+//----------------------------------------------------------------------
+int
+WaveformCableMappingTool::getChannelMapping(const EventContext& ctx, const Identifier id) const {
+  // Print where you are
+  ATH_MSG_DEBUG("in getChannelMapping()");
+  int channel = -1;
+
+  // Read Cond Handle
+  SG::ReadCondHandle<CondAttrListCollection> readHandle{m_readKey, ctx};
+  const CondAttrListCollection* readCdo{*readHandle}; 
+  if (readCdo==nullptr) {
+    ATH_MSG_FATAL("Null pointer to the read conditions object");
+    return channel;
+  }
+  // Get the validitiy range
+  EventIDRange rangeW;
+  if (not readHandle.range(rangeW)) {
+    ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key());
+    return channel;
+  }
+  ATH_MSG_DEBUG("Size of CondAttrListCollection " << readHandle.fullKey() << " readCdo->size()= " << readCdo->size());
+  ATH_MSG_DEBUG("Range of input is " << rangeW);
+  
+  // Read mapping info
+  CondAttrListCollection::const_iterator attrList{readCdo->begin()};
+  CondAttrListCollection::const_iterator end{readCdo->end()};
+  // CondAttrListCollection doesn't support C++11 type loops, no generic 'begin'
+  for (; attrList!=end; ++attrList) {
+    // A CondAttrListCollection is a map of ChanNum and AttributeList
+    CondAttrListCollection::ChanNum channelNumber{attrList->first};
+    const CondAttrListCollection::AttributeList &payload{attrList->second};
+    if (payload.exists("type") and not payload["type"].isNull()) {
+
+      std::string det_type{payload["type"].data<std::string>()};
+
+      int stationVal{payload["station"].data<int>()};
+      int plateVal  {payload["plate"].data<int>()};
+      int rowVal    {payload["row"].data<int>()};
+      int moduleVal {payload["module"].data<int>()};
+      int pmtVal    {payload["pmt"].data<int>()};
+      Identifier identifier;
+
+      // Ugh, cant use switch statement with strings
+      // Must do this using an if ladder
+      if (det_type == "calo") {
+	identifier = m_ecalID->pmt_id(rowVal, moduleVal, pmtVal);
+      }
+      else if (det_type == "veto") {
+	identifier = m_vetoID->pmt_id(stationVal, plateVal, pmtVal);
+      }
+      else if (det_type == "vetonu") {
+	identifier = m_vetoNuID->pmt_id(stationVal, plateVal, pmtVal);
+      }
+      else if (det_type == "trigger") {
+	identifier = m_triggerID->pmt_id(stationVal, plateVal, pmtVal);
+      }
+      else if (det_type == "preshower") {
+	identifier = m_preshowerID->pmt_id(stationVal, plateVal, pmtVal);
+      }
+      else if (det_type == "clock") {
+	// No valid identifiers for these
+	identifier = -1;
+      }
+      else if (det_type == "none") {
+	identifier = -1;
+      }
+      else {
+	ATH_MSG_WARNING("Detector type " << det_type << " not known for channel " << channelNumber << "!");
+	det_type = std::string("none");
+	identifier = -1;
+      }
+
+      // Is this the identifier we are looking for?
+      if (id != identifier) continue;
+
+      ATH_MSG_DEBUG("Mapped identifier " << det_type << " ID: " << identifier << " to digitizer channel " << channelNumber);
 
+      channel = channelNumber;
+      break;
+    }
+
+  } // End of loop over attributes
+
+  if (channel < 0)
+    ATH_MSG_WARNING("No channel found for identifier " << id << "!");
+
+  return channel;
+} 
+
+int
+WaveformCableMappingTool::getChannelMapping(const Identifier id) const {
+  const EventContext& ctx{Gaudi::Hive::currentContext()};
+  return getChannelMapping(ctx, id);
+}
 
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.h
index 19763d56e76bd1d5b18bcf74d20ad1b85d3c7acc..8443cb1815e4db7a9c0181515b5fd405cd27d48c 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.h
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformCableMappingTool.h
@@ -57,6 +57,11 @@ class WaveformCableMappingTool: public extends<AthAlgTool, IWaveformCableMapping
   virtual WaveformCableMap getCableMapping(const EventContext& ctx) const override;
   virtual WaveformCableMap getCableMapping(void) const override;
 
+  // Reverse mapping, reads idenfifier and returns digitizer channel
+  // Returns -1 if match not found for given identifier
+  virtual int getChannelMapping(const EventContext& ctx, const Identifier id) const override;
+  virtual int getChannelMapping(const Identifier id) const override;
+
  private:
   // Read Cond Handle
   SG::ReadCondHandleKey<CondAttrListCollection> m_readKey{this, "ReadKey", "/WAVE/DAQ/CableMapping", "Key of input cabling folder"};
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f163c65620a107ea8a5812818e5f8c5ba4049bac
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.cxx
@@ -0,0 +1,130 @@
+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS and FASER collaborations
+*/
+
+/** @file WaveformTimingTool.cxx Implementation file for WaveformTimingTool.
+    @author Eric Torrence (04/05/22)
+*/
+
+#include "WaveformTimingTool.h"
+
+//----------------------------------------------------------------------
+WaveformTimingTool::WaveformTimingTool (const std::string& type, const std::string& name, const IInterface* parent) :
+  base_class(type, name, parent)
+{
+}
+
+//----------------------------------------------------------------------
+StatusCode 
+WaveformTimingTool::initialize() {
+  // Read Cond Handle Key
+
+  ATH_MSG_DEBUG("WaveformTimingTool::initialize()");
+
+  ATH_CHECK(m_timingReadKey.initialize());
+  ATH_CHECK(m_offsetReadKey.initialize());
+
+  return StatusCode::SUCCESS;
+}
+
+//----------------------------------------------------------------------
+StatusCode
+WaveformTimingTool::finalize() {
+  // Print where you are
+  return StatusCode::SUCCESS;
+}
+
+//----------------------------------------------------------------------
+float 
+WaveformTimingTool::nominalTriggerTime(const EventContext& ctx) const {
+  // Print where you are
+  ATH_MSG_DEBUG("in nominalTriggerTime()");
+
+  float time=-1.;
+
+  // Read Cond Handle
+  SG::ReadCondHandle<AthenaAttributeList> readHandle{m_timingReadKey, ctx};
+  const AthenaAttributeList* readCdo(*readHandle);
+
+  if (readCdo==nullptr) {
+    ATH_MSG_FATAL("Null pointer to the read conditions object");
+    return time;
+  }
+
+  // Get the validitiy range
+  EventIDRange rangeW;
+  if (not readHandle.range(rangeW)) {
+    ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key());
+    return time;
+  }
+  ATH_MSG_DEBUG("Range of input is " << rangeW);
+  
+  // Read time info
+
+  const CondAttrListCollection::AttributeList &payload{*readCdo};
+  if (payload.exists("NominalTriggerTime") and not payload["NominalTriggerTime"].isNull()) {
+    time = payload["NominalTriggerTime"].data<float>();
+    ATH_MSG_DEBUG("Found nominal trigger time "<<time<<" ns");
+  } else {
+    ATH_MSG_WARNING("No valid nominal trigger time found!");
+  }
+
+  return time;
+} 
+
+//----------------------------------------------------------------------
+float
+WaveformTimingTool::triggerTimeOffset(const EventContext& ctx, int channel) const {
+
+  ATH_MSG_DEBUG("in triggerTimeOffset("<<channel<<")");
+
+  float time=0.;
+
+  // Read Cond Handle
+  SG::ReadCondHandle<CondAttrListCollection> readHandle{m_offsetReadKey, ctx};
+  const CondAttrListCollection* readCdo{*readHandle}; 
+  if (readCdo==nullptr) {
+    ATH_MSG_FATAL("Null pointer to the read conditions object");
+    return time;
+  }
+  // Get the validitiy range
+  EventIDRange rangeW;
+  if (not readHandle.range(rangeW)) {
+    ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key());
+    return time;
+  }
+  ATH_MSG_DEBUG("Size of CondAttrListCollection " << readHandle.fullKey() << " readCdo->size()= " << readCdo->size());
+  ATH_MSG_DEBUG("Range of input is " << rangeW);
+
+  // Read offset for specific channel
+  const CondAttrListCollection::AttributeList& payload{readCdo->attributeList(channel)};
+
+  if (payload.exists("TriggerOffset") and not payload["TriggerOffset"].isNull()) {
+    time = payload["TriggerOffset"].data<float>();
+    ATH_MSG_DEBUG("Found digitizer channel " << channel << " triger offset as " << time);
+  } else {
+    ATH_MSG_WARNING("No valid trigger offset found for channel "<<channel<<"!");
+  }
+
+  return time;
+
+}
+
+//----------------------------------------------------------------------
+float
+WaveformTimingTool::nominalTriggerTime(void) const {
+  const EventContext& ctx(Gaudi::Hive::currentContext());
+  return nominalTriggerTime(ctx);
+}
+
+float
+WaveformTimingTool::triggerTimeOffset(int channel) const {
+  const EventContext& ctx(Gaudi::Hive::currentContext());
+  return triggerTimeOffset(ctx, channel);
+}
+
+
+
+
+
+
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..b69e52b200907a369865f85b5b709aec5341d0dd
--- /dev/null
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformTimingTool.h
@@ -0,0 +1,69 @@
+// -*- C++ -*-
+
+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS and CERN collaborations
+*/
+
+/** @file WaveformTimingTool.h Header file for WaveformTimingTool.
+    @author Eric Torrence, 20/04/22
+*/
+
+// Multiple inclusion protection
+#ifndef WAVEFORM_TIMING_TOOL
+#define WAVEFORM_TIMING_TOOL
+
+// Include interface class
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "WaveformConditionsTools/IWaveformTimingTool.h"
+
+// Include Athena stuff
+#include "AthenaPoolUtilities/CondAttrListCollection.h"
+#include "StoreGate/ReadCondHandleKey.h"
+
+#include "GaudiKernel/ICondSvc.h"
+#include "Gaudi/Property.h"
+
+// Include Gaudi classes
+#include "GaudiKernel/EventContext.h"
+
+/** This class contains a Tool that reads Waveform timing data and makes it available 
+    to other algorithms. 
+*/
+
+class WaveformTimingTool: public extends<AthAlgTool, IWaveformTimingTool> {
+
+ public:
+  //----------Public Member Functions----------//
+  // Structors
+  WaveformTimingTool(const std::string& type, const std::string& name, const IInterface* parent); //!< Constructor
+  virtual ~WaveformTimingTool() = default; //!< Destructor
+
+  // Standard Gaudi functions
+  virtual StatusCode initialize() override; //!< Gaudi initialiser
+  virtual StatusCode finalize() override; //!< Gaudi finaliser
+
+  // Methods to return timing data
+  // Channels indexed by digitizer channel number
+  // All times are in ns
+
+  // Nominal trigger time (in ns) in the digitizer readout
+  virtual float nominalTriggerTime(void) const override;
+  virtual float nominalTriggerTime(const EventContext& ctx) const override;
+
+  // Channel-by-channel corrections to the nominal trigger time (in ns)
+  // A given channel should be centered at nominal + offset
+  virtual float triggerTimeOffset(int channel) const override;
+  virtual float triggerTimeOffset(const EventContext& ctx, int channel) const override;
+
+ private:
+
+  // Read Cond Handle
+  SG::ReadCondHandleKey<AthenaAttributeList> m_timingReadKey{this, "TimingReadKey", "/WAVE/DAQ/Timing", "Key of timing folder"};
+  SG::ReadCondHandleKey<CondAttrListCollection> m_offsetReadKey{this, "OffsetReadKey", "/WAVE/DAQ/TimingOffset", "Key of timing offset folder"};
+
+  ServiceHandle<ICondSvc> m_condSvc{this, "CondSvc", "CondSvc"};
+
+};
+
+//---------------------------------------------------------------------- 
+#endif // WAVEFORM_CABLE_MAPPING_TOOL
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/components/WaveformConditionsTools_entries.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/components/WaveformConditionsTools_entries.cxx
index 99371a2040253bcd84b16ff7a1886d922f2540af..f694fef647428802bfb98ab393e4f784180ae2b0 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/components/WaveformConditionsTools_entries.cxx
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/components/WaveformConditionsTools_entries.cxx
@@ -1,5 +1,7 @@
 #include "../WaveformRangeTool.h"
+#include "../WaveformTimingTool.h"
 #include "../WaveformCableMappingTool.h"
 
 DECLARE_COMPONENT( WaveformRangeTool )
+DECLARE_COMPONENT( WaveformTimingTool )
 DECLARE_COMPONENT( WaveformCableMappingTool )
diff --git a/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx
index df91c67cf77ef79b749934687ee910146ee3e309..59562e632c5704d53ca09eec67b0d036c6ac45b0 100644
--- a/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx
+++ b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx
@@ -27,6 +27,8 @@ namespace xAOD {
 
   AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, bcid_time, set_bcid_time )
 
+  AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, trigger_time, set_trigger_time )
+
   AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, raw_peak, set_raw_peak )
 
   AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, raw_integral, set_raw_integral )
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h
index 182003d90ab7f018f8b9afcccf45212b4cebc188..8e497633d1d673e6bb6220552874e5dbea1ce0e1 100644
--- a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h
@@ -74,6 +74,10 @@ namespace xAOD {
     float bcid_time() const;
     void set_bcid_time(float value);
 
+    /// Time with respect to nominal trigger time (including known offsets)
+    float trigger_time() const;
+    void set_trigger_time(float value);
+
     /// Raw values from waveform
     float raw_peak() const;
     void set_raw_peak(float value);