WaveformAnalysisLucyDDM.hh

#ifndef __RAT_WaveformAnalysisLucyDDM__
#define __RAT_WaveformAnalysisLucyDDM__
 
#include <TObject.h>
 
#include <RAT/DB.hh>
#include <RAT/DS/DigitPMT.hh>
#include <RAT/Digitizer.hh>
#include <RAT/FFTW1DTransformer.hh>
#include <RAT/Processor.hh>
#include <RAT/WaveformAnalyzerBase.hh>
#include <vector>
 
namespace RAT {
 
class WaveformAnalysisLucyDDM : public WaveformAnalyzerBase {
 public:
  WaveformAnalysisLucyDDM() : WaveformAnalysisLucyDDM("LucyDDM"){};
  WaveformAnalysisLucyDDM(std::string config_name) : WaveformAnalyzerBase("WaveformAnalysisLucyDDM", config_name) {
    Configure(config_name);
  };
  virtual ~WaveformAnalysisLucyDDM(){};
  void Configure(const std::string &config_name) override;
 
 protected:
  // Digitizer settings
  DBLinkPtr fDigit;
 
  // shaping parameters for SPE waveform
  double vpe_scale;   // Lognormal `m`
  double vpe_shape;   // Lognormal `sigma`
  double vpe_charge;  // nominal charge of a PE.
  double vpe_integral =
      -9999;  // to be overwritten by the charge * termohms. Set to invalid here because termOhms is not known yet.
  double ds;  // internal sampling period in ns.
  size_t vpe_nsamples;  // number of samples in the VPE waveform.
  std::vector<double> vpe_norm,
      vpe_norm_flipped;  // single PE waveform and its flipped version. Both normalized to integral 1.
  std::vector<std::complex<double>> vpe_norm_fft, vpe_norm_flipped_fft;  // and their fourier transforms.
 
  double roi_threshold;    // Waveform below this threshold will be zeroed out.
  double epsilon = 1e-10;  // Small value to avoid division by zero in deconvolution
 
  size_t max_iterations;     // max iterations to run in deconvolution
  double stopping_nll_diff;  // Stop if the change in the poisson negative log-likelihood is less than this value
 
  double peak_height_threshold;  // peak height threshold for finding hits
  double charge_threshold;       // remove a hit if it has charge below this threshold.
  double min_peak_distance;      // if peaks are closer than this, they will be merged.
 
  bool npe_estimate;  // if true, perform a final NPE estimation on all resolved wave packet in the deconvolved
                      // waveform. Estimate the number of PEs in each packet using a gaussian PDF on charge.
                      // The mean of the single PE charge distribution is centered around `vpe_charge`.
  bool npe_estimate_charge_width;  // the width of the gaussian single-PE charge distribution.
  size_t npe_estimate_max_pes;     // upper limit for the number of PEs in a single resolved wave packet.
 
  // FFT transfomration engines.
  std::unique_ptr<FFTW1DTransformer> fft = nullptr;
  std::unique_ptr<FFTW1DTransformer> ifft = nullptr;
 
  void DoAnalysis(DS::DigitPMT *pmt, const std::vector<UShort_t> &digitWfm) override;
 
  void RequestFFTSize(size_t size);
 
  std::vector<double> ConvolveFFT(const std::vector<double> &a, const std::vector<std::complex<double>> &b,
                                  size_t conv_size, double dt) const;
 
  std::vector<double> Deconvolve(const std::vector<double> &wfm, size_t &iterations_ran) const;
 
  void FindHits(const std::vector<double> &phi, std::vector<double> &out_times, std::vector<double> &out_charges,
                std::vector<double> &out_time_errors, std::vector<double> &out_charge_errors, double &chi2ndf) const;
 
  void MergeClosePeaks(std::vector<double> &times, std::vector<double> &charges, std::vector<double> &time_errors,
                       std::vector<double> &charge_errors);
 
  void ClampBelowThreshold(std::vector<double> &wfm, double thresh = -9999) const;
 
  std::vector<double> ReblurWaveform(const std::vector<double> &phi) const;
 
  std::vector<double> Resample(const std::vector<double> &wfm, size_t n_samples) const;
 
  double PoissonNLL(const std::vector<double> &wfm, const std::vector<double> &reblurred_wfm) const;
 
  double GaussianPulseTrain(double *x, double *p, size_t N) const;
};
 
}  // namespace RAT
#endif