Merge branch 'wavereco' into 'master'

Waveform reconstruction bugfix

See merge request faser/calypso!123

Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py

@@ -29,7 +29,7 @@ def WaveformClockRecCfg(flags, name="ClockRecAlg", **kwargs):
     acc = ComponentAccumulator()
 
     tool = ClockReconstructionTool(name="ClockReconstructionTool")
-
+    # tool.CheckResult = True
     kwargs.setdefault("ClockReconstructionTool", tool)
 
     recoAlg = CompFactory.ScintClockRecAlg(name, **kwargs)

Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx

@@ -64,7 +64,8 @@ ClockReconstructionTool::reconstruct(const ScintWaveform& raw_wave,
 
   ATH_MSG_DEBUG("Created double array with length " << wave.size() );
   ATH_MSG_DEBUG("First 10 elements:");
-  for (unsigned int i=0; i < std::min(10, N); i++)
+
+  for (int i=0; i < std::min(10, N); i++)
     ATH_MSG_DEBUG(" " << i << " " << wave[i]);
 
   // delta_nu = 1/T where T is the total waveform length
@@ -138,5 +139,34 @@ ClockReconstructionTool::reconstruct(const ScintWaveform& raw_wave,
 
   delete fftr2c;
 
+  // Check the result for sanity
+  if (clockdata->amplitude() < m_amplitude_min || clockdata->frequency() < m_frequency_min || clockdata->frequency() > m_frequency_max) {
+    ATH_MSG_WARNING("Clock reconstructed with frequency = " << clockdata->frequency() << " amplitude = " << clockdata->amplitude() << "!");
+  }
+
+  if (m_checkResult) checkResult(raw_wave, clockdata);
+
   return StatusCode::SUCCESS;
 }
+
+void
+ClockReconstructionTool::checkResult(const ScintWaveform& raw_wave,
+				     xAOD::WaveformClock* clockdata) const {
+
+  // Go through each element in raw_wave and make sure time in clockdata matches
+
+  float time;
+  for (unsigned int i=0; i<raw_wave.adc_counts().size(); i++) {
+    time = 2.*i; // Time in ns
+
+    float dt = clockdata->time_from_clock(time);
+
+    // Is raw_wave HI or LO?
+    bool hi = raw_wave.adc_counts()[i] > clockdata->dc_offset();
+
+    // Check for mismatch
+    if (((dt > 0.05) && (dt < 12.45) && !hi) || ((dt > 12.55) && (dt < 24.95) && hi) ) 
+      ATH_MSG_WARNING("Clock Time: " << time << " Time to edge: "  << dt << " found raw data: " << hi);
+  }
+
+}

Scintillator/ScintRecTools/src/ClockReconstructionTool.h

@@ -45,6 +45,17 @@ class ClockReconstructionTool: public extends<AthAlgTool, IClockReconstructionTo
   /// Minimum samples in the input waveform array to try FFT
   IntegerProperty m_minimumSamples{this, "MinimumSamples", 40};
 
+  /// Check reconstructed clock against waveform
+  BooleanProperty m_checkResult{this, "CheckResult", false};
+
+  void checkResult(const ScintWaveform& raw_wave,
+		   xAOD::WaveformClock* clockdata) const;
+
+  // Limits to print warnings
+  FloatProperty m_amplitude_min{this, "AmplitudeMin", 1000.};
+  FloatProperty m_frequency_min{this, "FrequencyMin", 40.0};
+  FloatProperty m_frequency_max{this, "FrequencyMax", 40.1};
+
 };
 
 #endif // SCINTRECTOOLS_CLOCKRECONSTRUCTIONTOOL_H

Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx

@@ -20,6 +20,7 @@
 
 #include <vector>
 #include <tuple>
+#include <math.h>
 
 // Constructor
 WaveformReconstructionTool::WaveformReconstructionTool(const std::string& type, const std::string& name, const IInterface* parent) :
@@ -168,8 +169,13 @@ WaveformReconstructionTool::reconstruct(const ScintWaveform& raw_wave,
     // Perform Gaussian fit to waveform
     WaveformFitResult gfit = fitGaussian(raw, wtime, wwave);
     if (! gfit.valid) {
-      ATH_MSG_WARNING( " Gaussian waveform fit failed! ");
-      hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED);
+      // 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);
+      }
     } 
 
     // Fit results (or raw if it failed)
@@ -455,6 +461,12 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons
 
   ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak);  
 
+  // Fix bad values
+  if (isnan(rfit.sigma)) {
+    ATH_MSG_DEBUG("Found sigma: " << rfit.sigma << " Replace with 2.");
+    rfit.sigma = 2.;
+  }
+
   return rfit;
 }
 
@@ -479,7 +491,7 @@ WaveformReconstructionTool::fitGaussian(const WaveformFitResult& raw, const std:
   gfunc.SetParameters(raw.peak, raw.mean, raw.sigma);
   gfunc.SetParError(0, std::sqrt(raw.peak));
   gfunc.SetParError(1, raw.sigma);
-  gfunc.SetParError(2, raw.sigma / 10.);
+  gfunc.SetParError(2, raw.sigma / 5.);
   gfunc.SetParLimits(2, 0., 20.);  // Constrain width
 
   TFitResultPtr gfitptr = tg.Fit(&gfunc, "QNS", "");
@@ -523,9 +535,9 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit,
   cbfit.peak = abs(gfit.peak);
   cbfit.mean = gfit.mean;
   cbfit.sigma = abs(gfit.sigma);
-  if (cbfit.sigma > 20.) cbfit.sigma = 5.;
-  cbfit.alpha = -0.25;  // Negative to get tail on high side
-  cbfit.nval = 3.;
+  if (cbfit.sigma > 20.) cbfit.sigma = 2.;
+  cbfit.alpha = -0.5;  // Negative to get tail on high side
+  cbfit.nval = 2.;
 
   // Define fit function and preset values
   TF1 cbfunc("cbfunc", "crystalball");
@@ -534,16 +546,27 @@ WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit,
   cbfunc.SetParameter(1, cbfit.mean); // Mean
   cbfunc.SetParError(1, cbfit.sigma);
   cbfunc.SetParameter(2, cbfit.sigma);  // Width
-  cbfunc.SetParError(2, cbfit.sigma/10.);
+  cbfunc.SetParError(2, cbfit.sigma/5.);
   cbfunc.SetParLimits(2, 0., 20.);
   cbfunc.SetParameter(3, cbfit.alpha); // Tail parameter (negative for high-side tail)
-  cbfunc.SetParError(3, -cbfit.alpha/10.);
+  cbfunc.SetParError(3, -cbfit.alpha/5.);
   cbfunc.SetParLimits(3, -10., 0.);
+  // Try fixing the tail first
+  cbfunc.FixParameter(4, cbfit.nval);
+
+  TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", "");
+
+  if (!cbfitptr->IsValid()) {
+    ATH_MSG_WARNING( " First Crystal Ball waveform fit failed! ");
+  }
+
+  // Now try releasing the tail parameter
+  cbfunc.ReleaseParameter(4);
   cbfunc.SetParameter(4, cbfit.nval);  // Tail power
-  cbfunc.SetParError(4, cbfit.nval/10.);
+  cbfunc.SetParError(4, cbfit.nval/5.);
   cbfunc.SetParLimits(4, 0., 1.E3);
 
-  TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", "");
+  cbfitptr = tg.Fit(&cbfunc, "QNS", "");
 
   cbfit.fit_status = cbfitptr;
   cbfit.valid = (cbfit.fit_status == 0);

Scintillator/ScintRecTools/src/WaveformReconstructionTool.h

@@ -80,7 +80,7 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
 
   //
   // Fraction of peak to set local hit time
-  FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.2};
+  FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.45};
 
   // Baseline algorithms
   WaveformBaselineData& findSimpleBaseline(const ScintWaveform& wave) const;

Tracker/TrackerRecAlgs/TrackerClusterFit/src/ClusterFitAlg.cxx

@@ -377,9 +377,9 @@ ClusterFitAlg::AddTrack(std::unique_ptr<TrackCollection>& tracks,
     s->push_back( GetState(fitResult, fitCovariance, nullptr, station) );
 
     for (const clusterInfo* cInfo : fitClusters)
-    {
-        s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
-    }
+      {
+	s->push_back( GetState(fitResult, fitCovariance, cInfo->cluster, station) );
+      }
 
     // Create and store track
     tracks->push_back(new Trk::Track(i, s , q));
@@ -495,4 +495,4 @@ ClusterFitAlg::Residuals(std::vector<const clusterInfo*>& fitClusters) const
     }
 }
 
-}
+}
\ No newline at end of file

xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h

@@ -105,7 +105,7 @@ namespace xAOD {
       this->set_status(this->status() | (1<<bit));
     }
     bool status_bit(WaveformStatus bit) const {
-      return (this->status() | (1<<bit));
+      return (this->status() & (1<<bit));
     }
 
     bool threshold() const {