util::SignalShaping Class Reference

#include <SignalShaping.h>

Public Member Functions
	SignalShaping ()
virtual	~SignalShaping ()
const std::vector< double > &	Response () const
const std::vector< double > &	Response_save () const
const std::vector< TComplex > &	ConvKernel () const
const std::vector< TComplex > &	Filter () const
const std::vector< TComplex > &	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< TComplex > &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< TComplex >	fConvKernel
std::vector< TComplex >	fFilter
std::vector< TComplex >	fDeconvKernel
int	fDeconvKernelPolarity
bool	fNorm

Detailed Description

Definition at line 69 of file SignalShaping.h.

Constructor & Destructor Documentation

util::SignalShaping::SignalShaping

(

)

Definition at line 19 of file SignalShaping.cxx.

  : fResponseLocked(false)
  , fFilterLocked  (false)
  , fNorm (true)
{}

util::SignalShaping::~SignalShaping ( )

virtual

Definition at line 29 of file SignalShaping.cxx.

{}

Member Function Documentation

void util::SignalShaping::AddFilterFunction

(

const std::vector< TComplex > &

filt

)

Definition at line 168 of file SignalShaping.cxx.

{
  // Make sure configuration is not locked.

  if(fFilterLocked)
    throw cet::exception("SignalShaping") << "Configuration locked.\n";

  // Get FFT service.

  art::ServiceHandle<util::LArFFT const> fft;

  // 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(fft->FFTSize() / 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::SignalShaping::AddResponseFunction	(	const std::vector< double > &	resp,
		bool	ResetResponse = false
	)

Definition at line 62 of file SignalShaping.cxx.

{
  // Make sure configuration is not locked.

  if(fResponseLocked)
    throw cet::exception("SignalShaping") << "Configuration locked.\n";

  // Get FFT service.

  art::ServiceHandle<util::LArFFT> fft;
  int nticks = fft->FFTSize();

  // 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);
    fft->DoFFT(fResponse, fConvKernel);
  }
  else {

    // Not the first response function.
    // Calculate the fourier transform of new response function.

    std::vector<TComplex> kern(nticks/2 + 1);
    fft->DoFFT(fResponse, kern);

    // Update overall convolution kernel.

    if (kern.size() != fConvKernel.size()) {
      throw cet::exception("SignalShaping") << __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.

    fft->DoInvFFT(fConvKernel, fResponse);
  }
}

void util::SignalShaping::CalculateDeconvKernel

(

)

const

Definition at line 251 of file SignalShaping.cxx.

{
  // Make sure configuration is not locked.

  if(fFilterLocked)
    throw cet::exception("SignalShaping") << "Configuration locked.\n";

  // Lock response configuration.

  LockResponse();

  // Get FFT service.

  art::ServiceHandle<util::LArFFT> fft;

  // Make sure filter function has been configured.

  if(fFilter.size() == 0)
    throw cet::exception("SignalShaping")
      << "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 = fft->FFTSize();
  if (2 * (fFilter.size() - 1) != n)
  if (fFilter.size() != fConvKernel.size()) {
    throw cet::exception("SignalShaping") << __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].Re()) <= 0.0001 && std::abs(fConvKernel[i].Im()) <= 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.);
  fft->DoInvFFT(const_cast<std::vector<TComplex>&>(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("SignalShaping") << __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("SignalShaping") << __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<TComplex>& util::SignalShaping::ConvKernel ( ) const

inline

Definition at line 79 of file SignalShaping.h.

{return fConvKernel;}

template<class T >

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

inline

Definition at line 167 of file SignalShaping.h.

{
  // Make sure response configuration is locked.
  if(!fResponseLocked)
    LockResponse();

  art::ServiceHandle<util::LArFFT> fft;

  // Make sure that time series has the correct size.
  if(int const n = func.size(); n != fft->FFTSize())
    throw cet::exception("SignalShaping") << "Bad time series size = " << n << "\n";

  fft->Convolute(func, const_cast<std::vector<TComplex>&>(fConvKernel));
}

const std::vector<TComplex>& util::SignalShaping::DeconvKernel ( ) const

inline

Definition at line 81 of file SignalShaping.h.

{return fDeconvKernel;}

template<class T >

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

inline

Definition at line 184 of file SignalShaping.h.

{
  // Make sure deconvolution kernel is configured.
  if(!fFilterLocked)
    CalculateDeconvKernel();

  art::ServiceHandle<util::LArFFT> fft;

  // Make sure that time series has the correct size.
  if(int const n = func.size(); n != fft->FFTSize())
    throw cet::exception("SignalShaping") << "Bad time series size = " << n << "\n";

  fft->Convolute(func, fDeconvKernel);
}

const std::vector<TComplex>& util::SignalShaping::Filter ( ) const

inline

Definition at line 80 of file SignalShaping.h.

{return fFilter;}

void util::SignalShaping::LockResponse

(

)

const

Definition at line 214 of file SignalShaping.cxx.

{
  // Do nothing if the response is already locked.

  if(!fResponseLocked) {

    // Get FFT service.

    art::ServiceHandle<util::LArFFT const> fft;

    // Make sure response has been configured.

    if(fResponse.size() == 0)
      throw cet::exception("SignalShaping")
        << "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 = fft->FFTSize();
    if (fResponse.size() != n)
      throw cet::exception("SignalShaping") << __func__ << ": inconsistent kernel size, "
        << fResponse.size() << " vs. " << n << "\n";
    if (2 * (fConvKernel.size() - 1) != n)
      throw cet::exception("SignalShaping") << __func__ << ": unexpected FFT size, "
        << n << " vs. expected " << (2 * (fConvKernel.size() - 1)) << "\n";

    // Set the lock flag.

    fResponseLocked = true;
  }
}

void util::SignalShaping::Reset

(

)

Definition at line 47 of file SignalShaping.cxx.

{
  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::SignalShaping::Response ( ) const

inline

Definition at line 77 of file SignalShaping.h.

{return fResponse;}

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

inline

Definition at line 78 of file SignalShaping.h.

{return fResponse_save;}

void util::SignalShaping::save_response ( )

inline

Definition at line 100 of file SignalShaping.h.

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

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

inline

Definition at line 101 of file SignalShaping.h.

{fNorm = flag;}

void util::SignalShaping::SetDeconvKernelPolarity

(

int

pol

)

Definition at line 197 of file SignalShaping.cxx.

{

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

  else
    fDeconvKernelPolarity = pol;

}

void util::SignalShaping::SetPeakResponseTime

(

double

tick

)

Definition at line 140 of file SignalShaping.cxx.

{
  // Make sure configuration is not locked.

  if(fResponseLocked)
    throw cet::exception("SignalShaping") << "Configuration locked.\n";

  // Get FFT service.

  art::ServiceHandle<util::LArFFT> fft;

  // Construct a delta-function response centered at tick zero.

  std::vector<double> delta(fft->FFTSize(), 0.);
  delta[0] = 1.;

  // Figure out peak of current overall response.

  double peak = fft->PeakCorrelation(delta, fResponse);

  // Shift peak response to desired tick.

  ShiftResponseTime(tick - peak);
}

void util::SignalShaping::ShiftResponseTime

(

double

ticks

)

Definition at line 117 of file SignalShaping.cxx.

{
  // Make sure configuration is not locked.

  if(fResponseLocked)
    throw cet::exception("SignalShaping") << "Configuration locked.\n";

  // Get FFT service.

  art::ServiceHandle<util::LArFFT> fft;

  // Update convolution kernel.

  fft->ShiftData(fConvKernel, ticks);

  // Recalculate overall response functiion.

  fft->DoInvFFT(fConvKernel, fResponse);
}

Member Data Documentation

std::vector<TComplex> util::SignalShaping::fConvKernel

private

Definition at line 146 of file SignalShaping.h.

std::vector<TComplex> util::SignalShaping::fDeconvKernel

mutableprivate

Definition at line 152 of file SignalShaping.h.

int util::SignalShaping::fDeconvKernelPolarity

private

Definition at line 157 of file SignalShaping.h.

std::vector<TComplex> util::SignalShaping::fFilter

private

Definition at line 149 of file SignalShaping.h.

bool util::SignalShaping::fFilterLocked

mutableprivate

Definition at line 139 of file SignalShaping.h.

bool util::SignalShaping::fNorm

private

Definition at line 160 of file SignalShaping.h.

std::vector<double> util::SignalShaping::fResponse

private

Definition at line 142 of file SignalShaping.h.

std::vector<double> util::SignalShaping::fResponse_save

private

Definition at line 143 of file SignalShaping.h.

bool util::SignalShaping::fResponseLocked

mutableprivate

Definition at line 138 of file SignalShaping.h.

---

The documentation for this class was generated from the following files:

• SignalShaping.h
• SignalShaping.cxx