#include "SignalShaper.h"

Public Member Functions
	SignalShaper (int fftsize, std::string fftopt)

virtual	~SignalShaper ()

const std::vector< double > &	Response () const

const std::vector< double > &	Response_save () const

const std::vector< std::complex< double > > &	ConvKernel () const

const std::vector< std::complex< double > > &	Filter () const

const std::vector< std::complex< double > > &	DeconvKernel () const

template<class T >
void	Convolute (std::vector< T > &func) const

template<class T >
void	Deconvolute (std::vector< T > &func) const

void	Reset ()

void	save_response ()

void	set_normflag (bool flag)

void	AddResponseFunction (const std::vector< double > &resp, bool ResetResponse=false)

void	ShiftResponseTime (double ticks)

void	SetPeakResponseTime (double tick)

void	AddFilterFunction (const std::vector< std::complex< double >> &filt)

void	SetDeconvKernelPolarity (int pol)

void	LockResponse () const

void	CalculateDeconvKernel () const

Private Attributes
bool	fResponseLocked

bool	fFilterLocked

std::vector< double >	fResponse

std::vector< double >	fResponse_save

std::vector< std::complex< double > >	fConvKernel

std::vector< std::complex< double > >	fFilter

std::vector< std::complex< double > >	fDeconvKernel

int	fDeconvKernelPolarity

bool	fNorm

int	fFFTSize

const void *	fPlan

const void *	rPlan

std::unique_ptr< util::LArFFTWPlan >	fFFTPlan

std::unique_ptr< util::LArFFTW >	fFFT

Detailed Description

Definition at line 13 of file SignalShaper.h.

Constructor & Destructor Documentation

util::SignalShaper::SignalShaper	(	int	fftsize,
		std::string	fftopt
	)

Definition at line 8 of file SignalShaper.cxx.

   : fResponseLocked(false)
   , fFilterLocked(false)
   , fNorm(true)
   , fFFTSize(fftsize)
   , fFFTPlan(new util::LArFFTWPlan(fftsize, fftopt))
   , fFFT(new util::LArFFTW(fftsize, fFFTPlan->fPlan, fFFTPlan->rPlan, 0))
 {}

util::SignalShaper::~SignalShaper ( )

virtual

Definition at line 20 of file SignalShaper.cxx.

20 {}

Member Function Documentation

void util::SignalShaper::AddFilterFunction ( const std::vector< std::complex< double >> & filt )

Definition at line 127 of file SignalShaper.cxx.

References fFFTSize, fFilter, fFilterLocked, filt, and n.

Referenced by set_normflag().

 {
   // Make sure configuration is not locked.
 
   if (fFilterLocked) throw cet::exception("SignalShaper") << "Configuration locked.\n";
 
   // If this is the first filter function, just copy the filter function.
   // Otherwise, update the overall filter function.
 
   if (fFilter.size() == 0) {
     fFilter = filt;
     fFilter.resize(fFFTSize / 2 + 1);
   }
   else {
     unsigned int n = std::min(fFilter.size(), filt.size());
     for (unsigned int i = 0; i < n; ++i)
       fFilter[i] *= filt[i];
     for (unsigned int i = n; i < fFilter.size(); ++i)
       fFilter[i] = 0.;
   }
 }

void util::SignalShaper::AddResponseFunction	(	const std::vector< double > &	resp,
		bool	ResetResponse = `false`
	)

Definition at line 38 of file SignalShaper.cxx.

References fConvKernel, fFFT, fFFTSize, fResponse, and fResponseLocked.

Referenced by set_normflag().

 {
   // Make sure configuration is not locked.
 
   if (fResponseLocked) throw cet::exception("SignalShaper") << "Configuration locked.\n";
 
   int nticks = fFFTSize;
 
   // Copy new response function into fResponse attribute, and pad or
   // truncate to correct size.
 
   fResponse = resp;
   fResponse.resize(nticks, 0.);
 
   // Is this the first response function?
 
   if (fConvKernel.size() == 0 || ResetResponse) {
 
     // This is the first response function.
     // Just calculate the fourier transform.
 
     fConvKernel.resize(nticks / 2 + 1);
     fFFT->DoFFT(fResponse, fConvKernel);
   }
   else {
 
     // Not the first response function.
     // Calculate the fourier transform of new response function.
 
     std::vector<std::complex<double>> kern(nticks / 2 + 1);
     fFFT->DoFFT(fResponse, kern);
 
     // Update overall convolution kernel.
 
     if (kern.size() != fConvKernel.size()) {
       throw cet::exception("SignalShaper") << __func__ << ": inconsistent kernel size, "
                                            << kern.size() << " vs. " << fConvKernel.size();
     }
     for (unsigned int i = 0; i < kern.size(); ++i)
       fConvKernel[i] *= kern[i];
 
     // Recalculate overall response function.
 
     fFFT->DoInvFFT(fConvKernel, fResponse);
   }
 }

void util::SignalShaper::CalculateDeconvKernel ( ) const

Definition at line 199 of file SignalShaper.cxx.

References util::abs(), fConvKernel, fDeconvKernel, fDeconvKernelPolarity, fFFT, fFFTSize, fFilter, fFilterLocked, fNorm, fResponse, LockResponse(), and n.

Referenced by set_normflag().

 {
   // Make sure configuration is not locked.
 
   if (fFilterLocked) throw cet::exception("SignalShaper") << "Configuration locked.\n";
 
   // Lock response configuration.
 
   LockResponse();
 
   // Make sure filter function has been configured.
 
   if (fFilter.size() == 0)
     throw cet::exception("SignalShaper") << "Filter function has not been configured.\n";
 
   // Make sure filter function has the correct size.
   // (Should always be the case if we get here.)
 
   unsigned int n = fFFTSize;
   if (2 * (fFilter.size() - 1) != n)
     if (fFilter.size() != fConvKernel.size()) {
       throw cet::exception("SignalShaper") << __func__ << ": inconsistent size, " << fFilter.size()
                                            << " vs. " << fConvKernel.size() << "\n";
     }
 
   // Calculate deconvolution kernel as the ratio of the
   // filter function and the convolution kernel.
 
   fDeconvKernel = fFilter;
   for (unsigned int i = 0; i < fDeconvKernel.size(); ++i) {
     if (std::abs(fConvKernel[i].real()) <= 0.0001 && std::abs(fConvKernel[i].imag()) <= 0.0001) {
       fDeconvKernel[i] = 0.;
     }
     else {
       fDeconvKernel[i] /= fConvKernel[i];
     }
   }
 
   // Normalize the deconvolution kernel.
 
   // Calculate the unnormalized deconvoluted response
   // (inverse FFT of filter function).
 
   std::vector<double> deconv(n, 0.);
   fFFT->DoInvFFT(const_cast<std::vector<std::complex<double>>&>(fFilter), deconv);
   //fFFT->DoInvFFT(fFilter, deconv);
 
   if (fNorm) {
     // Find the peak value of the response
     // Should normally be at zero, but don't assume that.
     // Use DeconvKernelPolairty to find what to normalize to
     double peak_response = 0;
     if (fDeconvKernelPolarity == -1) peak_response = 4096;
     for (unsigned int i = 0; i < fResponse.size(); ++i) {
       if ((fResponse[i] > peak_response) and (fDeconvKernelPolarity == 1))
         peak_response = fResponse[i];
       else if ((fResponse[i] < peak_response) and (fDeconvKernelPolarity == -1))
         peak_response = fResponse[i];
     }
     if (fDeconvKernelPolarity == -1) peak_response *= -1;
     if (peak_response <= 0.) {
       throw cet::exception("SignalShaper")
         << __func__ << ": peak should always be positive (got " << peak_response << ")\n";
     }
 
     // Find the peak value of the deconvoluted response
     // Should normally be at zero, but don't assume that.
 
     double peak_deconv = 0.;
     for (unsigned int i = 0; i < deconv.size(); ++i) {
       if (deconv[i] > peak_deconv) peak_deconv = deconv[i];
     }
     if (peak_deconv <= 0.) {
       throw cet::exception("SignalShaper")
         << __func__ << ": deconvolution peak should always be positive (got " << peak_deconv
         << ")\n";
     }
 
     // Multiply the deconvolution kernel by a factor such that
     // (Peak of response) = (Peak of deconvoluted response).
 
     double ratio = peak_response / peak_deconv;
     for (unsigned int i = 0; i < fDeconvKernel.size(); ++i)
       fDeconvKernel[i] *= ratio;
   }
 
   // Set the lock flag.
 
   fFilterLocked = true;
 }

const std::vector<std::complex<double> >& util::SignalShaper::ConvKernel ( ) const

inline

Definition at line 22 of file SignalShaper.h.

References fConvKernel.

22 { return fConvKernel; }

util::SignalShaper::fConvKernel

std::vector< std::complex< double > > fConvKernel

Definition: SignalShaper.h:87

template<class T >

void util::SignalShaper::Convolute ( std::vector< T > & func ) const

Referenced by DeconvKernel().

const std::vector<std::complex<double> >& util::SignalShaper::DeconvKernel ( ) const

inline

Definition at line 24 of file SignalShaper.h.

References Convolute(), Deconvolute(), fDeconvKernel, and Reset().

24 { return fDeconvKernel; }

util::SignalShaper::fDeconvKernel

std::vector< std::complex< double > > fDeconvKernel

Definition: SignalShaper.h:93

template<class T >

void util::SignalShaper::Deconvolute ( std::vector< T > & func ) const

Referenced by DeconvKernel().

const std::vector<std::complex<double> >& util::SignalShaper::Filter ( ) const

inline

Definition at line 23 of file SignalShaper.h.

References fFilter.

23 { return fFilter; }

util::SignalShaper::fFilter

std::vector< std::complex< double > > fFilter

Definition: SignalShaper.h:90

void util::SignalShaper::LockResponse ( ) const

Definition at line 167 of file SignalShaper.cxx.

References fConvKernel, fFFTSize, fResponse, fResponseLocked, and n.

Referenced by CalculateDeconvKernel(), and set_normflag().

 {
   // Do nothing if the response is already locked.
 
   if (!fResponseLocked) {
 
     // Make sure response has been configured.
 
     if (fResponse.size() == 0)
       throw cet::exception("SignalShaper") << "Response has not been configured.\n";
 
     // Make sure response and convolution kernel have the correct
     // size (should always be the case if we get here).
 
     unsigned int n = fFFTSize;
     if (fResponse.size() != n)
       throw cet::exception("SignalShaper")
         << __func__ << ": inconsistent kernel size, " << fResponse.size() << " vs. " << n << "\n";
     if (2 * (fConvKernel.size() - 1) != n)
       throw cet::exception("SignalShaper")
         << __func__ << ": unexpected FFT size, " << n << " vs. expected "
         << (2 * (fConvKernel.size() - 1)) << "\n";
 
     // Set the lock flag.
 
     fResponseLocked = true;
   }
 }

void util::SignalShaper::Reset ( )

Definition at line 24 of file SignalShaper.cxx.

References fConvKernel, fDeconvKernel, fDeconvKernelPolarity, fFilter, fFilterLocked, fResponse, and fResponseLocked.

Referenced by DeconvKernel().

 {
   fResponseLocked = false;
   fFilterLocked = false;
   fResponse.clear();
   fConvKernel.clear();
   fFilter.clear();
   fDeconvKernel.clear();
   //Set deconvolution polarity to + as default
   fDeconvKernelPolarity = +1;
 }

const std::vector<double>& util::SignalShaper::Response ( ) const

inline

Definition at line 20 of file SignalShaper.h.

References fResponse.

20 { return fResponse; }

util::SignalShaper::fResponse

std::vector< double > fResponse

Definition: SignalShaper.h:83

const std::vector<double>& util::SignalShaper::Response_save ( ) const

inline

Definition at line 21 of file SignalShaper.h.

References fResponse_save.

21 { return fResponse_save; }

util::SignalShaper::fResponse_save

std::vector< double > fResponse_save

Definition: SignalShaper.h:84

void util::SignalShaper::save_response ( )

inline

Definition at line 41 of file SignalShaper.h.

References fResponse, and fResponse_save.

     {
       fResponse_save.clear();
       fResponse_save = fResponse;
     }

void util::SignalShaper::set_normflag ( bool flag )

inline

Definition at line 46 of file SignalShaper.h.

References AddFilterFunction(), AddResponseFunction(), CalculateDeconvKernel(), fNorm, LockResponse(), SetDeconvKernelPolarity(), SetPeakResponseTime(), ShiftResponseTime(), and lar::dump::vector().

46 { fNorm = flag; }

util::SignalShaper::fNorm

bool fNorm

Definition: SignalShaper.h:101

void util::SignalShaper::SetDeconvKernelPolarity ( int pol )

Definition at line 151 of file SignalShaper.cxx.

References fDeconvKernelPolarity.

Referenced by set_normflag().

 {
 
   if ((pol != 1) and (pol != -1)) {
     throw cet::exception("SignalShaper")
       << __func__ << ": DeconvKernelPolarity should be +1 or -1 (got " << pol
       << "). Setting to +1\n";
     fDeconvKernelPolarity = +1;
   }
 
   else
     fDeconvKernelPolarity = pol;
 }

void util::SignalShaper::SetPeakResponseTime ( double tick )

Definition at line 105 of file SignalShaper.cxx.

References fFFT, fFFTSize, fResponse, fResponseLocked, and ShiftResponseTime().

Referenced by set_normflag().

 {
   // Make sure configuration is not locked.
 
   if (fResponseLocked) throw cet::exception("SignalShaper") << "Configuration locked.\n";
 
   // Construct a delta-function response centered at tick zero.
 
   std::vector<double> delta(fFFTSize, 0.);
   delta[0] = 1.;
 
   // Figure out peak of current overall response.
 
   double peak = fFFT->PeakCorrelation(delta, fResponse);
 
   // Shift peak response to desired tick.
 
   ShiftResponseTime(tick - peak);
 }

void util::SignalShaper::ShiftResponseTime ( double ticks )

Definition at line 88 of file SignalShaper.cxx.

References fConvKernel, fFFT, fResponse, and fResponseLocked.

Referenced by set_normflag(), and SetPeakResponseTime().

 {
   // Make sure configuration is not locked.
 
   if (fResponseLocked) throw cet::exception("SignalShaper") << "Configuration locked.\n";
 
   // Update convolution kernel.
 
   fFFT->ShiftData(fConvKernel, ticks);
 
   // Recalculate overall response functiion.
 
   fFFT->DoInvFFT(fConvKernel, fResponse);
 }