Add CB fit to waveform tool

Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.cxx

@@ -71,6 +71,9 @@ ScintWaveformRecAlg::execute(const EventContext& ctx) const {
     // Write something out
     ATH_MSG_INFO( "Channel " << wave->channel() << " baseline: " << *baseline );
 
+    // Reconstruct the hits
+    WaveformHitContainer* container = m_recoTool->reconstruct(*wave, *baseline);
+
     // This is transient, so delete it here
     delete baseline;
   }

Scintillator/ScintRecTools/ScintRecTools/IWaveformReconstructionTool.h

@@ -14,10 +14,14 @@
 
 // Base class
 #include "GaudiKernel/IAlgTool.h"
-#include "ScintRawEvent/ScintWaveform.h"
-#include "ScintRecEvent/WaveformBaselineData.h"
+//#include "ScintRawEvent/ScintWaveform.h"
+//#include "ScintRecEvent/WaveformBaselineData.h"
 #include "GaudiKernel/ToolHandle.h"
 
+class ScintWaveform;
+class WaveformBaselineData;
+class WaveformHitContainer;
+
 ///Interface for Waveform reco algorithms
 class IWaveformReconstructionTool : virtual public IAlgTool 
 {
@@ -28,12 +32,10 @@ class IWaveformReconstructionTool : virtual public IAlgTool
 
   virtual ~IWaveformReconstructionTool() = default;
 
+  // Find the baseline for a complete raw waveform
   virtual WaveformBaselineData* findBaseline(const ScintWaveform& wave) const = 0;
-
-  /** Reconstruct a single waveform into a list of hits
-   * @param[in] @c RDOs the raw data objects
-   */
-  //virtual WaveformHitCollection* reconstruct(const ScintWaveform& wave) const = 0;
+  // Reconstruct all peaks in a raw waveform
+  virtual WaveformHitContainer* reconstruct(const ScintWaveform& wave, const WaveformBaselineData& baseline) const = 0;
 
 };
 

Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx

@@ -14,6 +14,9 @@
 #include "TF1.h"
 #include "TFitResult.h"
 #include "TFitResultPtr.h"
+#include "TGraph.h"
+
+#include <vector>
 
 // Constructor
 WaveformReconstructionTool::WaveformReconstructionTool(const std::string& type, const std::string& name, const IInterface* parent) :
@@ -35,11 +38,167 @@ WaveformReconstructionTool::initialize() {
 }
 
 // Reconstruction step
-/*
-WaveformHitCollection*
-WaveformReconstructionTool::reconstruct(const ScintWaveform& wave) const {
+WaveformHitContainer*
+WaveformReconstructionTool::reconstruct(const ScintWaveform& rawWaveform, const WaveformBaselineData& baseline) const {
+
+  ATH_MSG_DEBUG(" reconstruct called ");
+
+  // Start out with empty WaveformHitContainer
+  WaveformHitContainer* container = new WaveformHitContainer();
+
+  // Create baseline-subtracted data array
+  std::vector<float> wave(rawWaveform.adc_counts().begin(), rawWaveform.adc_counts().end());
+  for (auto& element : wave)
+    element = baseline.mean - element;
+
+  // Define the fit functions here outside the loop
+  TF1 gfunc("gfunc", "gaus");
+  TF1 cbfunc("cbfunc", "crystalball");
+
+  // Now we iteratively find peaks and fit
+  while(true) {
+
+    // Find max value location in array
+    unsigned int imax = std::max_element(wave.begin(), wave.end()) - wave.begin();
+    float maxval = wave[imax];
+    ATH_MSG_DEBUG( "Found max value " << maxval << " at position " << imax );
+
+    // Check if this is over threshold
+    float pulseThreshold = 10.; // In baseline sigma
+    if (maxval < pulseThreshold*baseline.rms) {
+      ATH_MSG_DEBUG("Failed threshold");
+      break;
+    }
+
+    // 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 = (imax - 10) > 0 ? (imax - 10) : 0; 
+    int hi_edge = (imax + 40) < wave.size() ? (imax + 40) : wave.size();
+
+    ATH_MSG_DEBUG("Windowing waveform from " << lo_edge << " to " << hi_edge);
+    std::vector<float> window(wave.begin()+lo_edge, wave.begin()+hi_edge);
+
+    // Array of matching time values (in ns), find mean and rms of these data
+    std::vector<float> tvec(50, 0.);
+
+    double tot = 0.;
+    double sum = 0.;
+    double sum2 = 0.;
+    for (unsigned int i=0; i<tvec.size(); i++) {
+      tvec[i] = 2.*(lo_edge+i);
+      tot += window[i];
+      sum += tvec[i] * window[i];
+      sum2 += tvec[i] * tvec[i] * window[i];
+    }
+
+    // Initial parameters from window
+    double gmean = sum/tot;
+    double grms = std::sqrt(sum2/tot - gmean*gmean);
+    double gpeak = maxval;
+    ATH_MSG_DEBUG( "Initial Mean: " << gmean << " RMS: " << grms << " Peak: " << gpeak);  
+
+    // Start by creating a TGraph and fitting this with a Gaussian
+    // Must pass arrays to TGraph, use pointer to first vector element
+    TGraph tg(tvec.size(), &tvec[0], &window[0]);
+
+    // Define fit function and preset range
+    gfunc.SetParameters(gpeak, gmean, grms);
+    gfunc.SetParError(0, std::sqrt(gpeak));
+    gfunc.SetParError(1, grms);
+    gfunc.SetParError(2, 5.);
+    TFitResultPtr gfit = tg.Fit(&gfunc, "QNS", "");
+
+    //int gFitStatus = gfit;
+
+    if (!gfit->IsValid()) {
+      ATH_MSG_WARNING( " Gaussian waveform fit failed! ");
+      break;
+    } else {
+      // Improve estimation with fit results
+      gpeak = gfit->Parameter(0);
+      gmean = gfit->Parameter(1);
+      grms  = gfit->Parameter(2);
+      ATH_MSG_DEBUG("G Fit   Mean: " << gmean << " RMS: " << grms << " Peak: " << gpeak);
+    }
+
+    // Check for overflow
+    bool overflow = (rawWaveform.adc_counts()[imax] == 0);
+
+    if (overflow) {
+      ATH_MSG_INFO("Found waveform overflow");
+
+      // Want to remove overflowing points so fit works better
+      // Find first and last overflow here
+      // max_element will find the first overflow (imax), find the last
+      unsigned int iover;
+      for (iover = imax; iover < rawWaveform.adc_counts().size(); iover++) {
+	if (rawWaveform.adc_counts()[iover] != 0) break;
+      }
+
+      // Find the limits in the windowed waveform
+      unsigned int cut_lo = imax - lo_edge;
+      unsigned int cut_hi = iover - lo_edge;
+      if (cut_hi > window.size()) cut_hi = window.size();  // Don't cut beyond window
+
+      // Remove these entries
+      ATH_MSG_DEBUG("Removing entries from " << cut_lo << " to " << cut_hi);
+      tvec.erase(tvec.begin()+cut_lo, tvec.begin()+cut_hi);
+      window.erase(window.begin()+cut_lo, window.begin()+cut_hi);
+      ATH_MSG_DEBUG("Vector length now " << tvec.size());
+
+    }
+
+    // Now do crystal ball fit
+    double cbpeak = gpeak;
+    double cbmean = gmean;
+    double cbrms = grms;
+    if (cbrms < 0.) cbrms = abs(cbrms);
+    if (cbrms > 20.) cbrms = 5.;
+
+    double alpha = -0.25; // Negative puts tail on RH side
+    double nval = 3.;
+    // Preset everything to improve results
+    cbfunc.SetParameter(0, cbpeak); // Peak height
+    cbfunc.SetParError(0, std::sqrt(cbpeak));
+    cbfunc.SetParameter(1, cbmean); // Mean
+    cbfunc.SetParError(1, cbrms);
+    cbfunc.SetParameter(2, cbrms);  // Width
+    cbfunc.SetParError(2, 5.);
+    cbfunc.SetParLimits(2, 0., 20.);
+    cbfunc.SetParameter(3, alpha); // Tail parameter (negative for high-side tail)
+    cbfunc.SetParError(3, 0.05);
+    cbfunc.SetParLimits(3, -10., 0.);
+    cbfunc.SetParameter(4, nval);  // Tail power
+    cbfunc.SetParError(4, 1.);
+    cbfunc.SetParLimits(4, 0., 1.E3);
+    //cbfunc.SetParameters(peak, mean, rms, alpha, nval);
+
+    TFitResultPtr cbfit = tg.Fit(&cbfunc, "QNS", "");
+
+    double chi2 = cbfit->Chi2();
+    unsigned int ndf = cbfit->Ndf();
+    if (ndf == 0) ndf = 1;
+    double chi2ndf = chi2/ndf;
+
+    if (!cbfit->IsValid()) {
+      ATH_MSG_WARNING(" CrystalBall waveform fit failed!");
+      break;
+    } else {
+      // Improve estimation with fit results
+      cbpeak  = cbfit->Parameter(0);
+      cbmean  = cbfit->Parameter(1);
+      cbrms   = cbfit->Parameter(2);
+      alpha = cbfit->Parameter(3);
+      nval  = cbfit->Parameter(4);
+      ATH_MSG_DEBUG("CB Fit  Mean: " << cbmean << " RMS: " << cbrms << " Peak: " << cbpeak << " N: " << nval << " Alpha: " << alpha << " Chi2/Ndf: " << chi2ndf);
+    }
+  
+    break;
+  }
+
+  return container;
 }
-*/
+
 
 WaveformBaselineData*
 WaveformReconstructionTool::findBaseline(const ScintWaveform& wave) const {
@@ -55,7 +214,7 @@ WaveformBaselineData*
 WaveformReconstructionTool::findSimpleBaseline(const ScintWaveform& wave) const {
 
   ATH_MSG_DEBUG( "findSimpleBaseline called" );
-  ATH_MSG_DEBUG( wave );
+  //ATH_MSG_DEBUG( wave );
 
   // Calling algorithm checks for proper data
   // Here we just check algorithm-specific issues
@@ -91,11 +250,13 @@ WaveformReconstructionTool::findSimpleBaseline(const ScintWaveform& wave) const
 WaveformBaselineData*
 WaveformReconstructionTool::findAdvancedBaseline(const ScintWaveform& wave) const {
 
+  ATH_MSG_DEBUG( "findAdvancedBaseline called" );
+
   // First histogram to find most likely value
   // Turn these into configurables once this works
-  int nbins = 320;
+  int nbins = m_baselineRangeBins;
   double xlo = 0.;
-  double xhi = 16000.;
+  double xhi = m_baselineRange;
 
   TH1F h1("", "", nbins, xlo, xhi);
 
@@ -111,8 +272,8 @@ WaveformReconstructionTool::findAdvancedBaseline(const ScintWaveform& wave) cons
   ATH_MSG_DEBUG( "Found max bin at " << maxbinval << " counts");
 
   // Fill second histogram with integer resolution around the peak
-  nbins = 200;
-  xlo = int(maxbinval - 100.)-0.5;
+  nbins = m_baselineFitRange;
+  xlo = int(maxbinval - nbins/2)-0.5;
   xhi = xlo + nbins;
   ATH_MSG_DEBUG( "Filling 2nd histogram from " << xlo << " to " << xhi);
 
@@ -129,7 +290,7 @@ WaveformReconstructionTool::findAdvancedBaseline(const ScintWaveform& wave) cons
   ATH_MSG_DEBUG( "Initial Mean: " << mean << " RMS: " << rms << " Peak: " << peak );
 
   // Restrict range to +/- 2 sigma of mean
-  double window = 2.;  // Window range in sigma
+  double window = m_baselineFitWindow;  // Window range in sigma
   h2.GetXaxis()->SetRangeUser(mean-window*rms, mean+window*rms);
   mean = h2.GetMean();
   rms = h2.GetRMS();

Scintillator/ScintRecTools/src/WaveformReconstructionTool.h

@@ -12,7 +12,10 @@
 //Athena
 #include "AthenaBaseComps/AthAlgTool.h"
 #include "ScintRecTools/IWaveformReconstructionTool.h"
+
+#include "ScintRawEvent/ScintWaveform.h"
 #include "ScintRecEvent/WaveformBaselineData.h"
+#include "ScintRecEvent/WaveformHitContainer.h"
 
 //Gaudi
 #include "GaudiKernel/ToolHandle.h"
@@ -35,14 +38,27 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
   virtual WaveformBaselineData* findBaseline(const ScintWaveform& wave) const;
 
   /// Reconstruct hits from waveform
-  //virtual WaveformHitCollection* reconstruct(const ScintWaveform& wave) const;
+  virtual WaveformHitContainer* reconstruct(const ScintWaveform& wave, const WaveformBaselineData& baseline) const;
 
  private:
 
-  // Parameters
-  BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", true};
+  // Baseline Estimation Parameters
+  BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", false};
   IntegerProperty m_samplesForBaselineAverage{this, "SamplesForBaselineAverage", 40};
 
+  //
+  // Parameters of initial histogram to find baseline max
+  // Range and bins to use, ratio should be an integer (bin width)
+  IntegerProperty m_baselineRange{this, "BaselineRange", 16000};
+  IntegerProperty m_baselineRangeBins{this, "BaselineRangeBins", 320};
+
+  // 
+  // Parameters for the Gaussian fit to the baseline peak (in counts)
+  // Range is total range to use to find truncated mean and width
+  IntegerProperty m_baselineFitRange{this, "BaselineFitRange", 200};
+  // Fit window is value (in sigma) of trucated width to use in final fit
+  FloatProperty m_baselineFitWindow{this, "BaselineFitWindow", 2.};
+
   // Baseline algorithms
   WaveformBaselineData* findSimpleBaseline(const ScintWaveform& wave) const;
   WaveformBaselineData* findAdvancedBaseline(const ScintWaveform& wave) const;